{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#  第3章 k近邻法"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1．$k$近邻法是基本且简单的分类与回归方法。$k$近邻法的基本做法是：对给定的训练实例点和输入实例点，首先确定输入实例点的$k$个最近邻训练实例点，然后利用这$k$个训练实例点的类的多数来预测输入实例点的类。\n",
    "\n",
    "2．$k$近邻模型对应于基于训练数据集对特征空间的一个划分。$k$近邻法中，当训练集、距离度量、$k$值及分类决策规则确定后，其结果唯一确定。\n",
    "\n",
    "3．$k$近邻法三要素：距离度量、$k$值的选择和分类决策规则。常用的距离度量是欧氏距离及更一般的**pL**距离。$k$值小时，$k$近邻模型更复杂；$k$值大时，$k$近邻模型更简单。$k$值的选择反映了对近似误差与估计误差之间的权衡，通常由交叉验证选择最优的$k$。\n",
    "\n",
    "常用的分类决策规则是多数表决，对应于经验风险最小化。\n",
    "\n",
    "4．$k$近邻法的实现需要考虑如何快速搜索k个最近邻点。**kd**树是一种便于对k维空间中的数据进行快速检索的数据结构。kd树是二叉树，表示对$k$维空间的一个划分，其每个结点对应于$k$维空间划分中的一个超矩形区域。利用**kd**树可以省去对大部分数据点的搜索， 从而减少搜索的计算量。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 距离度量"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "设特征空间$x$是$n$维实数向量空间 ，$x_{i}, x_{j} \\in \\mathcal{X}$,$x_{i}=\\left(x_{i}^{(1)}, x_{i}^{(2)}, \\cdots, x_{i}^{(n)}\\right)^{\\mathrm{T}}$,$x_{j}=\\left(x_{j}^{(1)}, x_{j}^{(2)}, \\cdots, x_{j}^{(n)}\\right)^{\\mathrm{T}}$\n",
    "，则：$x_i$,$x_j$的$L_p$距离定义为:\n",
    "\n",
    "\n",
    "$L_{p}\\left(x_{i}, x_{j}\\right)=\\left(\\sum_{i=1}^{n}\\left|x_{i}^{(i)}-x_{j}^{(l)}\\right|^{p}\\right)^{\\frac{1}{p}}$\n",
    "\n",
    "- $p= 1$  曼哈顿距离\n",
    "- $p= 2$  欧氏距离\n",
    "- $p= \\infty$   切比雪夫距离"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import math\n",
    "from itertools import combinations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def L(x, y, p=2):\n",
    "    # x1 = [1, 1], x2 = [5,1]\n",
    "    if len(x) == len(y) and len(x) > 1:\n",
    "        sum = 0\n",
    "        for i in range(len(x)):\n",
    "            sum += math.pow(abs(x[i] - y[i]), p)\n",
    "        return math.pow(sum, 1 / p)\n",
    "    else:\n",
    "        return 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 课本例3.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "x1 = [1, 1]\n",
    "x2 = [5, 1]\n",
    "x3 = [4, 4]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(4.0, '1-[5, 1]')\n",
      "(4.0, '1-[5, 1]')\n",
      "(3.7797631496846193, '1-[4, 4]')\n",
      "(3.5676213450081633, '1-[4, 4]')\n"
     ]
    }
   ],
   "source": [
    "# x1, x2\n",
    "for i in range(1, 5):\n",
    "    r = {'1-{}'.format(c): L(x1, c, p=i) for c in [x2, x3]}\n",
    "    print(min(zip(r.values(), r.keys())))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "python实现，遍历所有数据点，找出$n$个距离最近的点的分类情况，少数服从多数"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "from sklearn.datasets import load_iris\n",
    "from sklearn.model_selection import train_test_split\n",
    "from collections import Counter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# data\n",
    "iris = load_iris()\n",
    "df = pd.DataFrame(iris.data, columns=iris.feature_names)\n",
    "df['label'] = iris.target\n",
    "df.columns = ['sepal length', 'sepal width', 'petal length', 'petal width', 'label']\n",
    "# data = np.array(df.iloc[:100, [0, 1, -1]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>sepal length</th>\n",
       "      <th>sepal width</th>\n",
       "      <th>petal length</th>\n",
       "      <th>petal width</th>\n",
       "      <th>label</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>5.1</td>\n",
       "      <td>3.5</td>\n",
       "      <td>1.4</td>\n",
       "      <td>0.2</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>4.9</td>\n",
       "      <td>3.0</td>\n",
       "      <td>1.4</td>\n",
       "      <td>0.2</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>4.7</td>\n",
       "      <td>3.2</td>\n",
       "      <td>1.3</td>\n",
       "      <td>0.2</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>4.6</td>\n",
       "      <td>3.1</td>\n",
       "      <td>1.5</td>\n",
       "      <td>0.2</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>5.0</td>\n",
       "      <td>3.6</td>\n",
       "      <td>1.4</td>\n",
       "      <td>0.2</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>145</th>\n",
       "      <td>6.7</td>\n",
       "      <td>3.0</td>\n",
       "      <td>5.2</td>\n",
       "      <td>2.3</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>146</th>\n",
       "      <td>6.3</td>\n",
       "      <td>2.5</td>\n",
       "      <td>5.0</td>\n",
       "      <td>1.9</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>147</th>\n",
       "      <td>6.5</td>\n",
       "      <td>3.0</td>\n",
       "      <td>5.2</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>148</th>\n",
       "      <td>6.2</td>\n",
       "      <td>3.4</td>\n",
       "      <td>5.4</td>\n",
       "      <td>2.3</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>149</th>\n",
       "      <td>5.9</td>\n",
       "      <td>3.0</td>\n",
       "      <td>5.1</td>\n",
       "      <td>1.8</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>150 rows × 5 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "     sepal length  sepal width  petal length  petal width  label\n",
       "0             5.1          3.5           1.4          0.2      0\n",
       "1             4.9          3.0           1.4          0.2      0\n",
       "2             4.7          3.2           1.3          0.2      0\n",
       "3             4.6          3.1           1.5          0.2      0\n",
       "4             5.0          3.6           1.4          0.2      0\n",
       "..            ...          ...           ...          ...    ...\n",
       "145           6.7          3.0           5.2          2.3      2\n",
       "146           6.3          2.5           5.0          1.9      2\n",
       "147           6.5          3.0           5.2          2.0      2\n",
       "148           6.2          3.4           5.4          2.3      2\n",
       "149           5.9          3.0           5.1          1.8      2\n",
       "\n",
       "[150 rows x 5 columns]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x24313189f48>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEHCAYAAACjh0HiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3df5wcdZ3n8df7hqyJCuSAuGsy0aB4eQghSyCCGBd/4C6K2RBREW7Z3Sgnq4sLHis+jOeioiconnjInS6gK4objFmMgAIqGAURuAnBRBLxJ25m4I4YNgHWoCF+7o+qSSZDz0zXdH+7q6rfz8djHtNdXV3zqWroT6rq8/l+FRGYmVnv+g/dDsDMzLrLicDMrMc5EZiZ9TgnAjOzHudEYGbW45wIzMx63D6p/4CkPmAAGIqIxaNeWwZcDAzliy6LiCvH295BBx0Uc+bMSRCpmVl9rV279tcRMaPRa8kTAXAOsAnYb4zXvxwR72h2Y3PmzGFgYKAtgZmZ9QpJvxrrtaSXhiT1A68Fxv1XvpmZdU/qewSfBN4N/H6cdV4vab2kVZJmN1pB0pmSBiQNbNmyJUmgZma9KlkikLQYeDgi1o6z2vXAnIiYD3wbuKrRShFxeUQsjIiFM2Y0vMRlZmaTlPIewSJgiaQTganAfpKujojTh1eIiK0j1r8C+GjCeMzMWrJz504GBwd54oknuh3KmKZOnUp/fz9Tpkxp+j3JEkFELAeWA0h6OfCukUkgX/7siHgof7qE7KaymVkpDQ4Osu+++zJnzhwkdTucp4gItm7dyuDgIAcffHDT7+t4H4GkCyQtyZ+eLek+ST8EzgaWdToeM7NmPfHEExx44IGlTAIAkjjwwAMLn7F0onyUiFgDrMkfnz9i+e6zBrO6Wb1uiItvvp8Ht+1g5vRpnHfCXJYumNXtsKxFZU0CwyYTX0cSgVmvWb1uiOXXbmDHzl0ADG3bwfJrNwA4GVjpeIgJswQuvvn+3Ulg2I6du7j45vu7FJHVxU033cTcuXM55JBDuOiii9qyTScCswQe3Laj0HKzZuzatYuzzjqLG2+8kY0bN7JixQo2btzY8nZ9acgsgZnTpzHU4Et/5vRpXYjGuqXd94nuvvtuDjnkEJ73vOcBcOqpp/K1r32NQw89tKU4fUZglsB5J8xl2pS+vZZNm9LHeSfM7VJE1mnD94mGtu0g2HOfaPW6oQnfO5ahoSFmz94zAEN/fz9DQ5Pf3jAnArMEli6YxYUnH86s6dMQMGv6NC48+XDfKO4hKe4TRcRTlrWjismXhswSWbpglr/4e1iK+0T9/f1s3rx59/PBwUFmzpw56e0N8xmBmVkCY90PauU+0Yte9CJ++tOf8stf/pLf/e53XHPNNSxZsmTiN07AicDMLIEU94n22WcfLrvsMk444QRe+MIXcsopp3DYYYe1GqovDZmZpTB8WbDd3eUnnngiJ554YjtC3M2JwMwskarcJ/KlITOzHudEYGbW45wIzMx6nBOBmVmPcyIwM+txTgTW81avG2LRRbdy8Hu+zqKLbm1pLBiz1N7ylrfwrGc9i3nz5rVtm04E1tNSDAxmltKyZcu46aab2rpNJwLraZ5AxpJavxIumQcfmJ79Xr+y5U0ed9xxHHDAAW0Ibg83lFlP8wQylsz6lXD92bAz/29p++bsOcD8U7oXVwM+I7CelmJgMDMAbrlgTxIYtnNHtrxknAisp3kCGUtm+2Cx5V3kS0PW01INDGbG/v3Z5aBGy0vGicB6XlUGBrOKOf78ve8RAEyZli1vwWmnncaaNWv49a9/TX9/Px/84Ac544wzWtqmE4F1Tbsn9jYrleEbwrdckF0O2r8/SwIt3ihesWJFG4LbmxOBdcVw/f5w6eZw/T7gZGD1Mf+U0lUINeKbxdYVrt83Kw8nAusK1+9bVUVEt0MY12TicyKwrnD9vlXR1KlT2bp1a2mTQUSwdetWpk6dWuh9vkdgXXHeCXP3ukcArt+38uvv72dwcJAtW7Z0O5QxTZ06lf7+YiWqTgTWFa7ftyqaMmUKBx98cLfDaDsnAusa1++blUPyRCCpDxgAhiJi8ajXngZ8ATgK2Aq8KSIeSB2TWdm4p8K6qRM3i88BNo3x2hnAv0XEIcAlwEc7EI9ZqXhOBOu2pIlAUj/wWuDKMVY5Cbgqf7wKOF6SUsZkVjbuqbBuS31G8Eng3cDvx3h9FrAZICKeBLYDB45eSdKZkgYkDZT5br3ZZLinwrotWSKQtBh4OCLWjrdag2VPKdCNiMsjYmFELJwxY0bbYjQrA/dUWLelPCNYBCyR9ABwDfBKSVePWmcQmA0gaR9gf+CRhDGZlY7nRLBuS5YIImJ5RPRHxBzgVODWiDh91GrXAX+dP35Dvk45W/bMElm6YBYXnnw4s6ZPQ8Cs6dO48OTDXTVkHdPxPgJJFwADEXEd8Fngi5J+RnYmcGqn4zErA/dUWDd1JBFExBpgTf74/BHLnwDe2IkYrHe8b/UGVty1mV0R9EmcdsxsPrz08G6HZVZa7iy2Wnnf6g1cfee/7n6+K2L3cycDs8Y8+qjVyoq7GswRO85yM3MisJrZNUatwVjLzcyJwGqmb4zG9LGWm5kTgdXMacfMLrTczHyz2Gpm+Iawq4bMmqeq9W8tXLgwBgYGuh2GmVmlSFobEQsbveYzAmurv7jiB3z/53tGCVn0/AP40luP7WJE3eM5BqwqfI/A2mZ0EgD4/s8f4S+u+EGXIuoezzFgVeJEYG0zOglMtLzOPMeAVYkTgVkCnmPAqsSJwCwBzzFgVeJEYG2z6PkHFFpeZ55jwKrEicDa5ktvPfYpX/q9WjXkOQasStxHYGbWA9xHYB2Tqna+yHZdv29WjBOBtc1w7fxw2eRw7TzQ0hdxke2misGsznyPwNomVe18ke26ft+sOCcCa5tUtfNFtuv6fbPinAisbVLVzhfZruv3zYpzIrC2SVU7X2S7rt83K843i61thm/Gtrtip8h2U8VgVmfuIzAz6wHuIyiRMtS4F42hDDGbWTpOBB1Uhhr3ojGUIWYzS8s3izuoDDXuRWMoQ8xmlpYTQQeVoca9aAxliNnM0nIi6KAy1LgXjaEMMZtZWk4EHVSGGveiMZQhZjNLyzeLO6gMNe5FYyhDzGaWlvsIzMx6QFf6CCRNBb4HPC3/O6si4v2j1lkGXAwM5Ysui4grU8Vkk/O+1RtYcddmdkXQJ3HaMbP58NLDW163LP0JZYnDrFtSXhr6LfDKiHhc0hTgdkk3RsSdo9b7ckS8I2Ec1oL3rd7A1Xf+6+7nuyJ2Px/9BV9k3bL0J5QlDrNumvBmsaSnSfrPkt4r6fzhn4neF5nH86dT8p9qXYcyVty1uenlRdYtS39CWeIw66Zmqoa+BpwEPAn8+4ifCUnqk3Qv8DDwrYi4q8Fqr5e0XtIqSbPH2M6ZkgYkDWzZsqWZP21tsmuMe0iNlhdZtyz9CWWJw6ybmrk01B8Rr57MxiNiF3CEpOnAVyXNi4gfjVjlemBFRPxW0tuAq4BXNtjO5cDlkN0snkwsNjl9UsMv8j6ppXVnTp/GUIMv2073J5QlDrNuauaM4A5Jje/2NSkitgFrgFePWr41In6bP70COKqVv2Ptd9oxDU/SGi4vsm5Z+hPKEodZN415RiBpA9k1/X2AN0v6BdkNYJHdApg/3oYlzQB2RsQ2SdOAVwEfHbXOsyPiofzpEmDTpPfEkhi+ydtMJVCRdcvSn1CWOMy6acw+AknPHe+NEfGrcTcszSe71NNHduaxMiIukHQBMBAR10m6kCwBPAk8Arw9In483nbdR2BmVtx4fQQTNpRJ+mJE/OVEyzql6okgVc16kfr9lNsusn9VPBaVs34l3HIBbB+E/fvh+PNh/indjsq6oNWGssNGbawPX8uflFQ160Xq91Nuu8j+VfFYVM76lXD92bAzvxm+fXP2HJwMbC9j3iyWtFzSY8B8SY/mP4+RlYJ+rWMR1kiqmvUi9fspt11k/6p4LCrnlgv2JIFhO3dky81GGDMRRMSFEbEvcHFE7Jf/7BsRB0bE8g7GWBupataL1O+n3HaR/avisaic7YPFllvPGu+M4EhJRwJfGX488qeDMdZGqrH9G9Xpj7c81baL7F8Vj0Xl7N9fbLn1rPH6CP5H/vO/gLvIGrquyB9fmj60+klVs16kfj/ltovsXxWPReUcfz5MGZVYp0zLlpuNMObN4oh4BYCka4AzI2JD/nwe8K7OhFcvqWrWi9Tvp9x2kf2r4rGonOEbwq4asgk0Uz56b0QcMdGyTql6+aiZWTe0Wj66SdKVwNVkncan4w7gnlKG3gCrOPczlFozieDNwNuBc/Ln3wM+nSwiK5Uy9AZYxbmfofQmHHQuIp6IiEsi4nX5zyUR8UQngrPuK0NvgFWc+xlKb7xB51ZGxCkjBp/by0SDzlk9lKE3wCrO/QylN96loeFLQYs7EYiVU5Hx+j22vzW0f392OajRciuF8TqLh4eHPh74g4j41cifzoRn3VaG3gCrOPczlF4zN4vnAKfnw1KvBW4DbouIe1MGZuVQht4Aqzj3M5TehH0Eu1fMJpd5K1kz2ayI6JvgLUm4j8DMrLiW+ggkvQ9YBDwTWEeWCG5ra4QllKoevsh2yzKuvnsDSqbuNfl1378iOnQsmrk0dDLZDGJfB74L3Fn38tFU9fBFtluWcfXdG1Ayda/Jr/v+FdHBY9FMH8GRZDeM7wb+FNgg6fa2RlEyqerhi2y3LOPquzegZOpek1/3/Suig8eimUtD84A/AV4GLAQ2U/NLQ6nq4Ytstyzj6rs3oGTqXpNf9/0rooPHYsIzAuCjwL5kQ0+/MCJeERG1rvtKNVZ+ke2WZVz9VMfCJqnucwzUff+K6OCxaObS0Gsj4mMRcUdE7Gx7BCWUqh6+yHbLMq6+ewNKpu41+XXfvyI6eCyauVncc1LVwxfZblnG1XdvQMnUvSa/7vtXRAePRdN9BGXhPgIzs+JanY/A2sj9CWYVccO5sPbzELtAfXDUMlj8ida3W8I+ifFGH72eBqOODouIJUkiqjH3J5hVxA3nwsBn9zyPXXuet5IMStonMealIUkvG++NEfHdJBFNoMqXhhZddGvD0TlnTZ/G99/zyo5s9/nLv9GwBLVP4ucXnjjpGMxq5YMHZF/+o6kP3v/I5Ld7ybwxRmKdDf/1R5PfbhMmdWmoW1/0deb+BLOKaJQExlverJL2SUxYPirpBZJWSdoo6RfDP50Irm7cn2BWERpjTM2xljerpH0SzTSU/RPZHMVPAq8AvgB8MWVQdeX+BLOKOGpZseXNKmmfRDOJYFpE3EJ2P+FXEfEBYPIXtHvY0gWzuPDkw5k1fRoiu4Z/4cmHt6U/odntfnjp4Zz+4ufsPgPokzj9xc/xjWKzkRZ/AhaesecMQH3Z81arhuafAn9+aXZPAGW///zSrlcNTdhHIOn7ZGMNrQJuBYaAiyKiK62lVb5ZbGbWLa32EbwTeDpwNvAhsrOBv27ij04Fvgc8Lf87qyLi/aPWeRrZpaajgK3AmyLigSZiKqxo/X7VxuAv0htQ92ORtE67SG15qjhS7l8Ja9zbpui+1flYjFJkhrL9gIiIx5pcX8AzIuJxSVOA24FzIuLOEev8LTA/It4m6VTgdRHxpvG2O5kzgtF19pBdQx/r8knR9bttdG/AsEaXfOp+LJ5Spw3ZNdh2nH6Pri0f1uiSQao4Uu5fym13W9F9q+GxGO+MoJmqoYWSNgDryeYi+KGkoyZ6X2Qez59OyX9GZ52TgKvyx6uA4/ME0lZFx9Sv2hj8ReYuqPuxSDqG+9rPN788VRwp96/OcwEU3bc6H4sGmrlZ/DngbyNiTkTMAc4iqySakKQ+SfcCDwPfioi7Rq0yi2x+AyLiSWA7cGCD7ZwpaUDSwJYtW5r503spWr9ftTH4i/QG1P1YJK3TLlJbniqOlPtX0hr3tii6b3U+Fg00kwgei4jdE9FExO1AU5eHImJXRBwB9ANH55PcjNToX/9P+faKiMsjYmFELJwxY0Yzf3ovRev3qzYGf5HegLofi6R12kVqy1PFkXL/Slrj3hZF963Ox6KBZhLB3ZL+UdLLJb1M0v8G1kg6UtKRzfyRiNgGrAFePeqlQWA2gKR9gP2BFvq3Gytav1+1MfiL9AbU/VgkrdMuUlueKo6U+1fSGve2KLpvdT4WDTRTNXRE/vv9o5a/hOxf7w17CiTNAHZGxDZJ04BXkc12NtJ1ZBVIPwDeANwaCcbFLjqmftXG4C8yd0Hdj0XSMdyHbwg3UzWUKo6U+1fnuQCK7ludj0UDyeYjkDSf7EZwH9mZx8qIuEDSBcBARFyXl5h+EVhAdiZwakSMO3yF+wjMzIprqY9A0h8CHwFmRsRrJB0KHBsRDero9oiI9WRf8KOXnz/i8RPAGyeKwczM0mnmHsHngZuBmfnzn5A1mdXa6nVDLLroVg5+z9dZdNGtrF431O2QrAzWr8yGEv7A9Oz3+pXtWTeVojGUYf+qtt0aaOYewUERsVLScsjKPCW1OBZruaWaQMYqrsikImWYgKRoDGXYv6pttyaaOSP4d0kHkpd1SnoxWb1/bVWuico6o0iTURkaklI2UVWtYa4Mn0eJNXNGcC5Zdc/z8wHoZpBV+NRW5ZqorDOKNBmVoSEpZRNV1RrmyvB5lNiEZwQRcQ/wMrJy0b8BDstvBNdW5ZqorDOKNBmVoSEpZRNV1RrmyvB5lFgzYw29kWxOgvuApcCXm20kq6rKNVFZZxRpMipDQ1LKJqqqNcyV4fMosWbuEfxDRDwm6aXACWS9AZ9OG1Z3pZpAxiquyKQiZZiApGgMZdi/qm23JpqZmGZdRCyQdCGwISL+eXhZZ0LcmxvKzMyKa3VimiFJ/0g+REQ+mUwzZxJmva3IJDZlUbWYyzJ5TFnimKRmEsEpZIPFfTwfN+jZwHlpwzKruNGT2MSuPc/L+sVatZjL0htQljha0EzV0G8i4tqI+Gn+/KGI+Gb60MwqrMgkNmVRtZjL0htQljha4Es8ZikUmcSmLKoWc1l6A8oSRwucCMxSKDKJTVlULeay9AaUJY4WOBGYpVBkEpuyqFrMZekNKEscLXAiMEth8Sdg4Rl7/jWtvux5GW+6DqtazGXpDShLHC1INjFNKu4jMDMrrtU+ArM0qlh7nSrmVPX7VTzG1nFOBNYdVay9ThVzqvr9Kh5j6wrfI7DuqGLtdaqYU9XvV/EYW1c4EVh3VLH2OlXMqer3q3iMrSucCKw7qlh7nSrmVPX7VTzG1hVOBNYdVay9ThVzqvr9Kh5j6wonAuuOKtZep4o5Vf1+FY+xdYX7CMzMesB4fQQ+IzBbvxIumQcfmJ79Xr+y89tNFYNZE9xHYL0tVa19ke263t+6zGcE1ttS1doX2a7r/a3LnAist6WqtS+yXdf7W5c5EVhvS1VrX2S7rve3LnMisN6Wqta+yHZd729d5kRgvS1VrX2R7bre37rMfQRmZj2gK30EkmZL+o6kTZLuk3ROg3VeLmm7pHvzH58LV10V6+Fd75+ej1uppewjeBL4+4i4R9K+wFpJ34qIjaPWuy0iFieMwzqlivXwrvdPz8et9JKdEUTEQxFxT/74MWATMCvV37MSqGI9vOv90/NxK72O3CyWNAdYANzV4OVjJf1Q0o2SDhvj/WdKGpA0sGXLloSRWkuqWA/vev/0fNxKL3kikPRM4F+Ad0bEo6Nevgd4bkT8MfApYHWjbUTE5RGxMCIWzpgxI23ANnlVrId3vX96Pm6llzQRSJpClgS+FBHXjn49Ih6NiMfzx98Apkg6KGVMllAV6+Fd75+ej1vppawaEvBZYFNENBxYXdIf5esh6eg8nq2pYrLEqlgP73r/9HzcSi9ZH4GklwK3ARuA3+eL3ws8ByAiPiPpHcDbySqMdgDnRsQd423XfQRmZsWN10eQrHw0Im4HNME6lwGXpYrBxrB+ZVaxsX0wu057/Pm9/a+zG86FtZ/PJotXXzZFZKuzg5lViOcj6DWu6d7bDefCwGf3PI9de547GViP8FhDvcY13Xtb+/liy81qyImg17ime2+xq9hysxpyIug1runem/qKLTerISeCXuOa7r0dtazYcrMaciLoNa7p3tviT8DCM/acAagve+4bxdZDPB+BmVkP6EofQS9ZvW6Ii2++nwe37WDm9Gmcd8Jcli6o0UCrde87qPv+lYGPcak5EbRo9bohll+7gR07syqToW07WH7tBoB6JIO69x3Uff/KwMe49HyPoEUX33z/7iQwbMfOXVx88/1diqjN6t53UPf9KwMf49JzImjRg9t2FFpeOXXvO6j7/pWBj3HpORG0aOb0aYWWV07d+w7qvn9l4GNcek4ELTrvhLlMm7J389G0KX2cd8LcLkXUZnXvO6j7/pWBj3Hp+WZxi4ZvCNe2amj4Zl5dKz7qvn9l4GNceu4jMDPrAeP1EfjSkFmdrV8Jl8yDD0zPfq9fWY1tW0f50pBZXaWs33dvQK34jMCsrlLW77s3oFacCMzqKmX9vnsDasWJwKyuUtbvuzegVpwIzOoqZf2+ewNqxYnArK5Szj3heS1qxX0EZmY9wH0EZmY2JicCM7Me50RgZtbjnAjMzHqcE4GZWY9zIjAz63FOBGZmPc6JwMysxyVLBJJmS/qOpE2S7pN0ToN1JOlSST+TtF7SkanisRZ43HmzWks5H8GTwN9HxD2S9gXWSvpWRGwcsc5rgBfkP8cAn85/W1l43Hmz2kt2RhARD0XEPfnjx4BNwOiJfE8CvhCZO4Hpkp6dKiabBI87b1Z7HblHIGkOsAC4a9RLs4DNI54P8tRkgaQzJQ1IGtiyZUuqMK0RjztvVnvJE4GkZwL/ArwzIh4d/XKDtzxlFLyIuDwiFkbEwhkzZqQI08bicefNai9pIpA0hSwJfCkirm2wyiAwe8TzfuDBlDFZQR533qz2UlYNCfgssCkiPjHGatcBf5VXD70Y2B4RD6WKySbB486b1V7KqqFFwF8CGyTdmy97L/AcgIj4DPAN4ETgZ8BvgDcnjMcma/4p/uI3q7FkiSAibqfxPYCR6wRwVqoYzMxsYu4sNjPrcU4EZmY9zonAzKzHORGYmfU4JwIzsx7nRGBm1uOcCMzMepyyUv7qkLQF+FW34xjDQcCvux1EQt6/6qrzvoH3rxnPjYiGg7VVLhGUmaSBiFjY7ThS8f5VV533Dbx/rfKlITOzHudEYGbW45wI2uvybgeQmPevuuq8b+D9a4nvEZiZ9TifEZiZ9TgnAjOzHudEMAmS+iStk3RDg9eWSdoi6d785790I8ZWSHpA0oY8/oEGr0vSpZJ+Jmm9pCO7EedkNLFvL5e0fcTnV6k5OSVNl7RK0o8lbZJ07KjXK/vZQVP7V9nPT9LcEXHfK+lRSe8ctU6Szy/lDGV1dg6wCdhvjNe/HBHv6GA8KbwiIsZqYHkN8IL85xjg0/nvqhhv3wBui4jFHYumvf4ncFNEvEHSHwBPH/V61T+7ifYPKvr5RcT9wBGQ/WMTGAK+Omq1JJ+fzwgKktQPvBa4stuxdNFJwBcicycwXdKzux1Ur5O0H3Ac2VzhRMTvImLbqNUq+9k1uX91cTzw84gYPYpCks/PiaC4TwLvBn4/zjqvz0/bVkma3aG42imAb0paK+nMBq/PAjaPeD6YL6uCifYN4FhJP5R0o6TDOhlci54HbAH+Kb90eaWkZ4xap8qfXTP7B9X9/EY6FVjRYHmSz8+JoABJi4GHI2LtOKtdD8yJiPnAt4GrOhJcey2KiCPJTkPPknTcqNcbzUVdlTrkifbtHrIxWf4Y+BSwutMBtmAf4Ejg0xGxAPh34D2j1qnyZ9fM/lX58wMgv+S1BPhKo5cbLGv583MiKGYRsETSA8A1wCslXT1yhYjYGhG/zZ9eARzV2RBbFxEP5r8fJrtGefSoVQaBkWc6/cCDnYmuNRPtW0Q8GhGP54+/AUyRdFDHA52cQWAwIu7Kn68i++IcvU4lPzua2L+Kf37DXgPcExH/r8FrST4/J4ICImJ5RPRHxByyU7dbI+L0keuMul63hOymcmVIeoakfYcfA38G/GjUatcBf5VXMLwY2B4RD3U41MKa2TdJfyRJ+eOjyf4f2drpWCcjIv4vsFnS3HzR8cDGUatV8rOD5vavyp/fCKfR+LIQJPr8XDXUBpIuAAYi4jrgbElLgCeBR4Bl3YxtEv4Q+Gr+/9I+wD9HxE2S3gYQEZ8BvgGcCPwM+A3w5i7FWlQz+/YG4O2SngR2AKdGtdrv/w74Un554RfAm2vy2Q2baP8q/flJejrwp8DfjFiW/PPzEBNmZj3Ol4bMzHqcE4GZWY9zIjAz63FOBGZmPc6JwMysxzkRmBWUj3DZaOTZhsvb8PeWSjp0xPM1kmo7Ubt1nhOBWfktBQ6dcC2zSXIisNrJO4i/ng889iNJb8qXHyXpu/mAczcPd4Hn/8L+pKQ78vWPzpcfnS9bl/+eO97fbRDD5yT9n/z9J+XLl0m6VtJNkn4q6WMj3nOGpJ/k8Vwh6TJJLyHrUL9Y2Rj1z89Xf6Oku/P1/6RNh856lDuLrY5eDTwYEa8FkLS/pClkg5CdFBFb8uTw34G35O95RkS8JB+E7nPAPODHwHER8aSkVwEfAV7fZAz/jWwIkrdImg7cLenb+WtHAAuA3wL3S/oUsAv4B7Kxcx4DbgV+GBF3SLoOuCEiVuX7A7BPRBwt6UTg/cCrJnOgzMCJwOppA/BxSR8l+wK9TdI8si/3b+VfpH3AyDFaVgBExPck7Zd/ee8LXCXpBWQjPE4pEMOfkQ1Q+K78+VTgOfnjWyJiO4CkjcBzgYOA70bEI/nyrwD/aZztX5v/XgvMKRCX2VM4EVjtRMRPJB1FNibLhZK+STbS6H0RcexYb2vw/EPAdyLidZLmAGsKhCHg9fmsU3sWSseQnQkM20X2/2Gj4YXHM7yN4febTZrvEVjtSJoJ/CYirgY+Tna55X5ghvI5biVN0d6TlgzfR3gp2YiO24H9yaYLhOKDB94M/N2IkTAXTLD+3cDLJP1HSfuw9yWox8jOTsyS8L8krI4OJ7u5+ntgJ/D2iPidpDcAlz3SjnEAAACjSURBVEran+y//U8C9+Xv+TdJd5DNQz183+BjZJeGziW7Zl/Eh/Ltr8+TwQPAmPPoRsSQpI8Ad5GNL78R2J6/fA1whaSzyUbXNGsrjz5qPU/SGuBdETHQ5TieGRGP52cEXwU+FxGjJy83aztfGjIrjw9IupdsspxfUsFpFq2afEZgZtbjfEZgZtbjnAjMzHqcE4GZWY9zIjAz63FOBGZmPe7/A5LawLRyzHcuAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(df[:50]['sepal length'], df[:50]['sepal width'], label='0')\n",
    "plt.scatter(df[50:100]['sepal length'], df[50:100]['sepal width'], label='1')\n",
    "plt.xlabel('sepal length')\n",
    "plt.ylabel('sepal width')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = np.array(df.iloc[:100, [0, 1, -1]])\n",
    "X, y = data[:,:-1], data[:,-1]\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "class KNN:\n",
    "    def __init__(self, X_train, y_train, n_neighbors=3, p=2):\n",
    "        \"\"\"\n",
    "        parameter: n_neighbors 临近点个数\n",
    "        parameter: p 距离度量\n",
    "        \"\"\"\n",
    "        self.n = n_neighbors\n",
    "        self.p = p\n",
    "        self.X_train = X_train\n",
    "        self.y_train = y_train\n",
    "\n",
    "    def predict(self, X):\n",
    "        # 取出n个点\n",
    "        knn_list = []\n",
    "        for i in range(self.n):\n",
    "            dist = np.linalg.norm(X - self.X_train[i], ord=self.p)\n",
    "            knn_list.append((dist, self.y_train[i]))\n",
    "\n",
    "        for i in range(self.n, len(self.X_train)):\n",
    "            max_index = knn_list.index(max(knn_list, key=lambda x: x[0]))\n",
    "            dist = np.linalg.norm(X - self.X_train[i], ord=self.p)\n",
    "            if knn_list[max_index][0] > dist:\n",
    "                knn_list[max_index] = (dist, self.y_train[i])\n",
    "\n",
    "        # 统计\n",
    "        knn = [k[-1] for k in knn_list]\n",
    "        count_pairs = Counter(knn)\n",
    "#         max_count = sorted(count_pairs, key=lambda x: x)[-1]\n",
    "        max_count = sorted(count_pairs.items(), key=lambda x: x[1])[-1][0]\n",
    "        return max_count\n",
    "\n",
    "    def score(self, X_test, y_test):\n",
    "        right_count = 0\n",
    "        n = 10\n",
    "        for X, y in zip(X_test, y_test):\n",
    "            label = self.predict(X)\n",
    "            if label == y:\n",
    "                right_count += 1\n",
    "        return right_count / len(X_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "clf = KNN(X_train, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "clf.score(X_test, y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Test Point: 1.0\n"
     ]
    }
   ],
   "source": [
    "test_point = [6.0, 3.0]\n",
    "print('Test Point: {}'.format(clf.predict(test_point)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x2431324bf88>"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEHCAYAAACjh0HiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3df5xVdb3v8dfHcY6AIhTSOcKg44/yYQLJD3+FZqYd/EH+IDVNTlHeOHrsaJmUZMcflKnhTa95s4NmmSJKHMIflVoqqZl4BxhBJRNLZQbvlShQBD2An/vHWgPDds/MXrP3d++11n4/H4957Nlrr73ms9bW/WGt9fl8v+buiIhI/dqh1gGIiEhtKRGIiNQ5JQIRkTqnRCAiUueUCERE6pwSgYhIndsx9B8wswagBWh39wkFr00GZgDt8aIb3f2W7ra32267eXNzc4BIRUTya9GiRX9198HFXgueCIALgOXArl28fre7f7nUjTU3N9PS0lKRwERE6oWZvdLVa0EvDZlZE3AC0O2/8kVEpHZC3yO4Hvg68G4363zazJaa2VwzG1ZsBTObYmYtZtayevXqIIGKiNSrYInAzCYAr7v7om5Wuw9odveRwG+B24qt5O4z3X2su48dPLjoJS4REemlkPcIxgEnmtnxQB9gVzO7w90ndazg7ms6rX8zcE3AeEQkhTZt2kRbWxtvv/12rUPJhT59+tDU1ERjY2PJ7wmWCNx9GjANwMw+DlzUOQnEy3d399fipycS3VQWkTrS1tZG//79aW5uxsxqHU6muTtr1qyhra2Nvfbaq+T3Vb2PwMymm9mJ8dPzzew5M3sGOB+YXO14RKS23n77bQYNGqQkUAFmxqBBgxKfXVWjfBR3XwAsiH+/tNPyrWcNInkzf0k7Mx58gVVrNzJkYF+mjt+Pk0cNrXVYqaQkUDm9OZZVSQQi9Wb+knamzVvGxk1bAGhfu5Fp85YBKBlI6miICZEAZjz4wtYk0GHjpi3MePCFGkUk0jUlApEAVq3dmGi5lG7WLGhuhh12iB5nzSpve2vXruWHP/xhr957/fXXs2HDhvICKHDppZfy29/+ttt1FixYwJNPPlmxv6lEIBLAkIF9Ey2X0syaBVOmwCuvgHv0OGVKeckgbYlg+vTpHHPMMd2uo0QgkgFTx+9H38aG7Zb1bWxg6vj9ahRRPlxyCRR+727YEC3vrYsvvpiXXnqJAw88kKlTpzJjxgwOOuggRo4cyWWXXQbAW2+9xQknnMBHPvIRhg8fzt13380NN9zAqlWrOOqoozjqqKO63P4uu+zC1772NUaPHs3RRx9Nx+gIra2tHHrooYwcOZJTTjmFv//97wBMnjyZuXPnAtHYapdddhmjR49mxIgR/PGPf+Tll1/mRz/6Eddddx0HHnggjz/+eO93PqZEIBLAyaOGctXEEQwd2BcDhg7sy1UTR+hGcZlefTXZ8lJcffXV7LPPPrS2tvLJT36SF198kaeffprW1lYWLVrEY489xgMPPMCQIUN45plnePbZZzn22GM5//zzGTJkCI8++iiPPvpol9t/6623GD16NIsXL+bII4/kiiuuAOBzn/sc11xzDUuXLmXEiBFblxfabbfdWLx4Meeeey7XXnstzc3NnHPOOXz1q1+ltbWVI444ovc7H1PVkEggJ48aqi/+Cttjj+hyULHllfDQQw/x0EMPMWrUKADWr1/Piy++yBFHHMFFF13EN77xDSZMmJDoy3eHHXbgM5/5DACTJk1i4sSJrFu3jrVr13LkkUcC8PnPf57TTjut6PsnTpwIwJgxY5g3b145u9d1jEG2KiISwJVXQr9+2y/r1y9aXgnuzrRp02htbaW1tZUVK1Zw9tln86EPfYhFixYxYsQIpk2bxvTp03v9N5LW+e+0004ANDQ0sHnz5l7/3e4oEYhIZpx1FsycCXvuCWbR48yZ0fLe6t+/P2+++SYA48eP59Zbb2X9+vUAtLe38/rrr7Nq1Sr69evHpEmTuOiii1i8ePF73tuVd999d+s1/zvvvJPDDz+cAQMG8L73vW/r9f3bb79969lB0pgrQZeGRCRTzjqrvC/+QoMGDWLcuHEMHz6c4447js9+9rMcdthhQHSj94477mDFihVMnTqVHXbYgcbGRm666SYApkyZwnHHHcfuu+/e5X2CnXfemeeee44xY8YwYMAA7r77bgBuu+02zjnnHDZs2MDee+/NT37yk5Jj/tSnPsWpp57KPffcww9+8IOy7xOYu5e1gWobO3asa4YykfxYvnw5+++/f63DCGaXXXbZeoZRLcWOqZktcvexxdbXpSERkTqnS0MiIhVwyCGH8M4772y37Pbbb6/62UBvKBGIiFTAwoULax1Cr+nSkIhInVMiEBGpc7o0JHVPE8hIvdMZgdS1jglk2tduxNk2gcz8Je21Dk2q6IEHHmC//fZj33335eqrr651OFWnRCB1TRPIyJYtWzjvvPP49a9/zfPPP8/s2bN5/vnnax1WVenSkNQ1TSCTPZW+lPf000+z7777svfeewNwxhlncM899/DhD3+4UiGnns4IpK5pAplsCXEpr729nWHDhm193tTURHt7fV0aVCKQuqYJZLIlxKW8YsPsJB0hNOt0aUjqWsclBVUNZUOIS3lNTU2sXLly6/O2tjaGDBnS6+1lkRKB1D1NIJMdQwb2pb3Il345l/IOOuggXnzxRf7yl78wdOhQ7rrrLu68885ywswcXRqSmpm/pJ1xVz/CXhf/knFXP6KSTelRiEt5O+64IzfeeCPjx49n//335/TTT+eAAw4oN9RM0RmB1ETHTb+O670dN/0A/etcuhTqUt7xxx/P8ccfX4kQM0mJQGqiu5t+SgTSHV3KqzxdGpKaUP2+SHooEUhNqH5fJD2UCKQmVL8vkh66RyA1ofp9kfRQIpCa0U0/kXQIfmnIzBrMbImZ3V/ktZ3M7G4zW2FmC82sOXQ8Immknora+uIXv8gHPvABhg8fXutQaqIa9wguAJZ38drZwN/dfV/gOuCaKsQjkiqaE6H2Jk+ezAMPPFDrMGomaCIwsybgBOCWLlY5Cbgt/n0ucLTV22hPUvc0J0JCS+fAdcPh8oHR49I5ZW/yYx/7GO9///srEFw2hT4juB74OvBuF68PBVYCuPtmYB0wqHAlM5tiZi1m1rJ69epQsYrUhHoqElg6B+47H9atBDx6vO/8iiSDehYsEZjZBOB1d1/U3WpFlr1nTFh3n+nuY9197ODBgysWo0gaqKcigYenw6aCBLlpY7Rcei3kGcE44EQzexm4C/iEmd1RsE4bMAzAzHYEBgB/CxiTSOqopyKBdW3JlktJgiUCd5/m7k3u3gycATzi7pMKVrsX+Hz8+6nxOu+dJUIkx04eNZSrJo5g6MC+GDB0YF+umjhCpbXFDGhKtlxKUvU+AjObDrS4+73Aj4HbzWwF0ZnAGdWORyQN1FNRoqMvje4JdL481Ng3Wl6GM888kwULFvDXv/6VpqYmrrjiCs4+++wyg82OqiQCd18ALIh/v7TT8reB06oRg9SPb81fxuyFK9niToMZZx4yjO+cPKLWYUkljDw9enx4enQ5aEBTlAQ6lvfS7NmzKxBcdqmzWHLlW/OXccdTr259vsV963Mlg5wYeXrZX/yyPQ06J7kye+HKRMtFRIlAcmZLF7UGXS2XdFCNSOX05lgqEUiuNHTRmN7Vcqm9Pn36sGbNGiWDCnB31qxZQ58+fRK9T/cIJFfOPGTYdvcIOi+XdGpqaqKtrQ2NGlAZffr0oakpWTmtEoHkSscNYVUNZUdjYyN77bVXrcOoa5a107GxY8d6S0tLrcMQEckUM1vk7mOLvaYzAqmos27+A79/adsoIeP2eT+zvnRYDSOqnflL2jUDm2SCbhZLxRQmAYDfv/Q3zrr5DzWKqHY0x4BkiRKBVExhEuhpeZ5pjgHJEiUCkQA0x4BkiRKBSACaY0CyRIlAKmbcPsWn+utqeZ5pjgHJEiUCqZhZXzrsPV/69Vo1pDkGJEvURyAiUgfURyBVE6p2Psl2Vb8vkowSgVRMR+18R9lkR+08UNYXcZLthopBJM90j0AqJlTtfJLtqn5fJDklAqmYULXzSbar+n2R5JQIpGJC1c4n2a7q90WSUyKQiglVO59ku6rfF0lON4ulYjpuxla6YifJdkPFIJJn6iMQEakD6iNIkTTUuCeNIQ0xi0g4SgRVlIYa96QxpCFmEQlLN4urKA017kljSEPMIhKWEkEVpaHGPWkMaYhZRMJSIqiiNNS4J40hDTGLSFhKBFWUhhr3pDGkIWYRCUs3i6soDTXuSWNIQ8wiEpb6CERE6kBN+gjMrA/wGLBT/HfmuvtlBetMBmYA7fGiG939llAxSe98a/4yZi9cyRZ3Gsw485BhfOfkEWWvm5b+hLTEIVIrIS8NvQN8wt3Xm1kj8ISZ/drdnypY7253/3LAOKQM35q/jDueenXr8y3uW58XfsEnWTct/QlpiUOklnq8WWxmO5nZZ83sm2Z2acdPT+/zyPr4aWP8k63rUMLshStLXp5k3bT0J6QlDpFaKqVq6B7gJGAz8Fannx6ZWYOZtQKvA79x94VFVvu0mS01s7lmNqyL7UwxsxYza1m9enUpf1oqZEsX95CKLU+yblr6E9ISh0gtlXJpqMndj+3Nxt19C3CgmQ0EfmFmw9392U6r3AfMdvd3zOwc4DbgE0W2MxOYCdHN4t7EIr3TYFb0i7zBrKx1hwzsS3uRL9tq9yekJQ6RWirljOBJMyt+t69E7r4WWAAcW7B8jbu/Ez+9GRhTzt+RyjvzkKInaUWXJ1k3Lf0JaYlDpJa6PCMws2VE1/R3BL5gZn8mugFsRLcARna3YTMbDGxy97Vm1hc4BrimYJ3d3f21+OmJwPJe74kE0XGTt5RKoCTrpqU/IS1xiNRSl30EZrZnd29091e63bDZSKJLPQ1EZx5z3H26mU0HWtz9XjO7iigBbAb+Bpzr7n/sbrvqIxARSa67PoIeG8rM7HZ3/5eellVL1hNBqJr1JPX7IbedZP+yeCwyZ+kceHg6rGuDAU1w9KUw8vRaRyU1UG5D2QEFG2tA1/J7JVTNepL6/ZDbTrJ/WTwWmbN0Dtx3PmyKb4avWxk9ByUD2U6XN4vNbJqZvQmMNLM34p83iUpB76lahDkSqmY9Sf1+yG0n2b8sHovMeXj6tiTQYdPGaLlIJ10mAne/yt37AzPcfdf4p7+7D3L3aVWMMTdC1awnqd8Pue0k+5fFY5E569qSLZe61d0ZwWgzGw38vOP3zj9VjDE3Qo3tX6xOv7vlobadZP+yeCwyZ0BTsuVSt7rrI/if8c//BhYSNXTdHP9+Q/jQ8idUzXqS+v2Q206yf1k8Fplz9KXQWJBYG/tGy0U66fJmsbsfBWBmdwFT3H1Z/Hw4cFF1wsuXUDXrSer3Q247yf5l8VhkTscNYVUNSQ9KKR9tdfcDe1pWLVkvHxURqYVyy0eXm9ktwB1EncaTUAdwXUlDb4BknPoZUq2URPAF4Fzggvj5Y8BNwSKSVElDb4BknPoZUq/HQefc/W13v87dT4l/rnP3t6sRnNReGnoDJOPUz5B63Q06N8fdT+80+Nx2ehp0TvIhDb0BknHqZ0i97i4NdVwKmlCNQCSdkozXr7H9pagBTdHloGLLJRW66yzuGB76aOAf3P2Vzj/VCU9qLQ29AZJx6mdIvVJuFjcDk+JhqRcBjwOPu3tryMAkHdLQGyAZp36G1Ouxj2DritHkMl8iaiYb6u4NPbwlCPURiIgkV1YfgZl9CxgH7AIsIUoEj1c0whQKVQ+fZLtpGVdfvQEpk/ea/LzvXxJVOhalXBqaSDSD2C+B3wFP5b18NFQ9fJLtpmVcffUGpEzea/Lzvn9JVPFYlNJHMJrohvHTwCeBZWb2REWjSJlQ9fBJtpuWcfXVG5Ayea/Jz/v+JVHFY1HKpaHhwBHAkcBYYCU5vzQUqh4+yXbTMq6+egNSJu81+XnfvySqeCx6PCMArgH6Ew09vb+7H+Xuua77CjVWfpLtpmVc/VDHQnop73MM5H3/kqjisSjl0tAJ7v49d3/S3TdVPIIUClUPn2S7aRlXX70BKZP3mvy8718SVTwWpdwsrjuh6uGTbDct4+qrNyBl8l6Tn/f9S6KKx6LkPoK0UB+BiEhy5c5HIBWk/gSRjLj/Qlj0U/AtYA0wZjJM+H75201hn0R3o4/eR5FRRzu4+4lBIsox9SeIZMT9F0LLj7c99y3bnpeTDFLaJ9HlpSEzO7K7N7r774JE1IMsXxoad/UjRUfnHDqwL7+/+BNV2e4+035VtAS1wYyXrjq+1zGI5MoV74++/AtZA1z2t95v97rhXYzEOgy++mzvt1uCXl0aqtUXfZ6pP0EkI4olge6WlyqlfRI9lo+a2QfNbK6ZPW9mf+74qUZweaP+BJGMsC7G1OxqealS2idRSkPZT4jmKN4MHAX8DLg9ZFB5pf4EkYwYMznZ8lKltE+ilETQ190fJrqf8Iq7Xw70/oJ2HTt51FCumjiCoQP7YkTX8K+aOKIi/Qmlbvc7J49g0qF7bD0DaDBj0qF76EaxSGcTvg9jz952BmAN0fNyq4ZGng6fuiG6J4BFj5+6oeZVQz32EZjZ74nGGpoLPAK0A1e7e01aS7N8s1hEpFbK7SP4CtAPOB/4NtHZwOdL+KN9gMeAneK/M9fdLytYZyeiS01jgDXAZ9z95RJiSixp/X7WxuBP0huQ92MRtE47SW15qDgSbHfWLLjkEnj1VdhjD7jySjjrrMpsO3OS7luej0WBJDOU7Qq4u79Z4voG7Ozu682sEXgCuMDdn+q0zr8BI939HDM7AzjF3T/T3XZ7c0ZQWGcP0TX0ri6fJF2/1gp7AzoUu+ST92PxnjptiK7BVuL0u7C2vEOxSwah4kiw3VmzYMoU2LBh27J+/WDmzC6SQchjV2tJ9y2Hx6K7M4JSqobGmtkyYCnRXATPmNmYnt7nkfXx08b4pzDrnATcFv8+Fzg6TiAVlXRM/ayNwZ9k7oK8H4ugY7gv+mnpy0PFkWC7l1yyfRKA6Pkll5S/7cxJum95PhZFlHKz+Fbg39y92d2bgfOIKol6ZGYNZtYKvA78xt0XFqwylGh+A9x9M7AOGFRkO1PMrMXMWlavXl3Kn95O0vr9rI3Bn6Q3IO/HImiddpLa8lBxJNjuq+89Sex2eVpr3Csi6b7l+VgUUUoieNPdt05E4+5PACVdHnL3Le5+INAEHBxPctNZsX/9v+fby91nuvtYdx87ePDgUv70dpLW72dtDP4kvQF5PxZB67ST1JaHiiPBdvfYo/iqXS1Pa417RSTdtzwfiyJKSQRPm9l/mtnHzexIM/shsMDMRpvZ6FL+iLuvBRYAxxa81AYMAzCzHYEBQBn928Ulrd/P2hj8SXoD8n4sgtZpJ6ktDxVHgu1eeWV0T6Czfv2i5eVuO3OS7luej0URpVQNHRg/Xlaw/KNE/3ov2lNgZoOBTe6+1sz6AscQzXbW2b1EFUh/AE4FHvEA42InHVM/a2PwJ5m7IO/HIugY7h03hEupGgoVR4LtdtwQLrlqKM9zASTdtzwfiyKCzUdgZiOJbgQ3EJ15zHH36WY2HWhx93vjEtPbgVFEZwJnuHu3w1eoj0BEJLmy+gjM7B+B7wJD3P04M/swcJi7F6mj28bdlxJ9wRcuv7TT728Dp/UUg4iIhFPKPYKfAg8CQ+LnfyJqMsu1+UvaGXf1I+x18S8Zd/UjzF/SXuuQJA2WzomGEr58YPS4dE5l1g0laQxp2L+sbTcHSrlHsJu7zzGzaRCVeZpZmWOxpluoCWQk45JMKpKGCUiSxpCG/cvadnOilDOCt8xsEHFZp5kdSlTvn1uZa6KS6kjSZJSGhqSQTVQpaJhLxXZzopQzgguJqnv2iQegG0xU4ZNbmWuikupI0mSUhoakkE1UKWiYS8V2c6LHMwJ3XwwcSVQu+q/AAfGN4NzKXBOVVEeSJqM0NCSFbKJKQcNcKrabE6WMNXQa0ZwEzwEnA3eX2kiWVZlropLqSNJklIaGpJBNVClomEvFdnOilHsE/+Hub5rZ4cB4ot6Am8KGVVuhJpCRjEsyqUgaJiBJGkMa9i9r282JUiamWeLuo8zsKmCZu9/Zsaw6IW5PDWUiIsmVOzFNu5n9J/EQEfFkMqWcSYjUtyST2KRF1mJOy+QxaYmjl0pJBKcTDRZ3bTxu0O7A1LBhiWRc4SQ2vmXb87R+sWYt5rT0BqQljjKUUjW0wd3nufuL8fPX3P2h8KGJZFiSSWzSImsxp6U3IC1xlEGXeERCSDKJTVpkLea09AakJY4yKBGIhJBkEpu0yFrMaekNSEscZVAiEAkhySQ2aZG1mNPSG5CWOMqgRCASwoTvw9izt/1r2hqi52m86dohazGnpTcgLXGUIdjENKGoj0BEJLly+whEwshi7XWomEPV72fxGEvVKRFIbWSx9jpUzKHq97N4jKUmdI9AaiOLtdehYg5Vv5/FYyw1oUQgtZHF2utQMYeq38/iMZaaUCKQ2shi7XWomEPV72fxGEtNKBFIbWSx9jpUzKHq97N4jKUmlAikNrJYex0q5lD1+1k8xlIT6iMQEakD3fUR6IxAZOkcuG44XD4welw6p/rbDRWDSAnURyD1LVStfZLtqt5fakxnBFLfQtXaJ9mu6v2lxpQIpL6FqrVPsl3V+0uNKRFIfQtVa59ku6r3lxpTIpD6FqrWPsl2Ve8vNaZEIPUtVK19ku2q3l9qTH0EIiJ1oCZ9BGY2zMweNbPlZvacmV1QZJ2Pm9k6M2uNf3QunHVZrIdXvX94Om6pFrKPYDPwNXdfbGb9gUVm9ht3f75gvcfdfULAOKRaslgPr3r/8HTcUi/YGYG7v+bui+Pf3wSWA0ND/T1JgSzWw6vePzwdt9Srys1iM2sGRgELi7x8mJk9Y2a/NrMDunj/FDNrMbOW1atXB4xUypLFenjV+4en45Z6wROBme0C/BfwFXd/o+DlxcCe7v4R4AfA/GLbcPeZ7j7W3ccOHjw4bMDSe1msh1e9f3g6bqkXNBGYWSNREpjl7vMKX3f3N9x9ffz7r4BGM9stZEwSUBbr4VXvH56OW+qFrBoy4MfAcncvOrC6mf1TvB5mdnAcz5pQMUlgWayHV71/eDpuqResj8DMDgceB5YB78aLvwnsAeDuPzKzLwPnElUYbQQudPcnu9uu+ghERJLrro8gWPmouz8BWA/r3AjcGCoG6cLSOVHFxrq26Drt0ZfW97/O7r8QFv00mizeGqIpIsudHUwkQzQfQb1RTff27r8QWn687blv2fZcyUDqhMYaqjeq6d7eop8mWy6SQ0oE9UY13dvzLcmWi+SQEkG9UU339qwh2XKRHFIiqDeq6d7emMnJlovkkBJBvVFN9/YmfB/Gnr3tDMAaoue6USx1RPMRiIjUgZr0EdST+UvamfHgC6xau5EhA/sydfx+nDwqRwOt5r3vIO/7lwY6xqmmRFCm+UvamTZvGRs3RVUm7Ws3Mm3eMoB8JIO89x3kff/SQMc49XSPoEwzHnxhaxLosHHTFmY8+EKNIqqwvPcd5H3/0kDHOPWUCMq0au3GRMszJ+99B3nfvzTQMU49JYIyDRnYN9HyzMl730He9y8NdIxTT4mgTFPH70ffxu2bj/o2NjB1/H41iqjC8t53kPf9SwMd49TTzeIyddwQzm3VUMfNvLxWfOR9/9JAxzj11EcgIlIHuusj0KUhkTxbOgeuGw6XD4wel87JxralqnRpSCSvQtbvqzcgV3RGIJJXIev31RuQK0oEInkVsn5fvQG5okQgklch6/fVG5ArSgQieRWyfl+9AbmiRCCSVyHnntC8FrmiPgIRkTqgPgIREemSEoGISJ1TIhARqXNKBCIidU6JQESkzikRiIjUOSUCEZE6p0QgIlLngiUCMxtmZo+a2XIze87MLiiyjpnZDWa2wsyWmtnoUPFIGTTuvEiuhZyPYDPwNXdfbGb9gUVm9ht3f77TOscBH4x/DgFuih8lLTTuvEjuBTsjcPfX3H1x/PubwHKgcCLfk4CfeeQpYKCZ7R4qJukFjTsvkntVuUdgZs3AKGBhwUtDgZWdnrfx3mSBmU0xsxYza1m9enWoMKUYjTsvknvBE4GZ7QL8F/AVd3+j8OUib3nPKHjuPtPdx7r72MGDB4cIU7qicedFci9oIjCzRqIkMMvd5xVZpQ0Y1ul5E7AqZEySkMadF8m9kFVDBvwYWO7u3+9itXuBz8XVQ4cC69z9tVAxSS9o3HmR3AtZNTQO+BdgmZm1xsu+CewB4O4/An4FHA+sADYAXwgYj/TWyNP1xS+SY8ESgbs/QfF7AJ3XceC8UDGIiEjP1FksIlLnlAhEROqcEoGISJ1TIhARqXNKBCIidU6JQESkzikRiIjUOYtK+bPDzFYDr9Q6ji7sBvy11kEEpP3LrjzvG2j/SrGnuxcdrC1ziSDNzKzF3cfWOo5QtH/Zled9A+1fuXRpSESkzikRiIjUOSWCyppZ6wAC0/5lV573DbR/ZdE9AhGROqczAhGROqdEICJS55QIesHMGsxsiZndX+S1yWa22sxa45//UYsYy2FmL5vZsjj+liKvm5ndYGYrzGypmY2uRZy9UcK+fdzM1nX6/DI1J6eZDTSzuWb2RzNbbmaHFbye2c8OStq/zH5+ZrZfp7hbzewNM/tKwTpBPr+QM5Tl2QXAcmDXLl6/292/XMV4QjjK3btqYDkO+GD8cwhwU/yYFd3tG8Dj7j6hatFU1v8CHnD3U83sH4B+Ba9n/bPraf8go5+fu78AHAjRPzaBduAXBasF+fx0RpCQmTUBJwC31DqWGjoJ+JlHngIGmtnutQ6q3pnZrsDHiOYKx93/293XFqyW2c+uxP3Li6OBl9y9cBSFIJ+fEkFy1wNfB97tZp1Px6dtc81sWJXiqiQHHjKzRWY2pcjrQ4GVnZ63xcuyoKd9AzjMzJ4xs1+b2QHVDK5MewOrgZ/Ely5vMbOdC9bJ8mdXyv5Bdj+/zs4AZhdZHuTzUyJIwMwmAK+7+6JuVrsPaHb3kcBvgduqElxljXP30USnoRWHANkAAASASURBVOeZ2ccKXi82F3VW6pB72rfFRGOyfAT4ATC/2gGWYUdgNHCTu48C3gIuLlgny59dKfuX5c8PgPiS14nAz4u9XGRZ2Z+fEkEy44ATzexl4C7gE2Z2R+cV3H2Nu78TP70ZGFPdEMvn7qvix9eJrlEeXLBKG9D5TKcJWFWd6MrT0765+xvuvj7+/VdAo5ntVvVAe6cNaHP3hfHzuURfnIXrZPKzo4T9y/jn1+E4YLG7/78irwX5/JQIEnD3ae7e5O7NRKduj7j7pM7rFFyvO5HopnJmmNnOZta/43fgn4FnC1a7F/hcXMFwKLDO3V+rcqiJlbJvZvZPZmbx7wcT/T+yptqx9oa7/19gpZntFy86Gni+YLVMfnZQ2v5l+fPr5EyKXxaCQJ+fqoYqwMymAy3ufi9wvpmdCGwG/gZMrmVsvfCPwC/i/5d2BO509wfM7BwAd/8R8CvgeGAFsAH4Qo1iTaqUfTsVONfMNgMbgTM8W+33/w7Mii8v/Bn4Qk4+uw497V+mPz8z6wd8EvjXTsuCf34aYkJEpM7p0pCISJ1TIhARqXNKBCIidU6JQESkzikRiIjUOSUCkYTiES6LjTxbdHkF/t7JZvbhTs8XmFluJ2qX6lMiEEm/k4EP97iWSC8pEUjuxB3Ev4wHHnvWzD4TLx9jZr+LB5x7sKMLPP4X9vVm9mS8/sHx8oPjZUvix/26+7tFYrjVzP5P/P6T4uWTzWyemT1gZi+a2fc6vedsM/tTHM/NZnajmX2UqEN9hkVj1O8Tr36amT0dr39EhQ6d1Cl1FkseHQuscvcTAMxsgJk1Eg1CdpK7r46Tw5XAF+P37OzuH40HobsVGA78EfiYu282s2OA7wKfLjGGS4iGIPmimQ0Enjaz38avHQiMAt4BXjCzHwBbgP8gGjvnTeAR4Bl3f9LM7gXud/e58f4A7OjuB5vZ8cBlwDG9OVAioEQg+bQMuNbMriH6An3czIYTfbn/Jv4ibQA6j9EyG8DdHzOzXeMv7/7AbWb2QaIRHhsTxPDPRAMUXhQ/7wPsEf/+sLuvAzCz54E9gd2A37n73+LlPwc+1M3258WPi4DmBHGJvIcSgeSOu//JzMYQjclylZk9RDTS6HPuflhXbyvy/NvAo+5+ipk1AwsShGHAp+NZp7YtNDuE6Eygwxai/w+LDS/cnY5tdLxfpNd0j0Byx8yGABvc/Q7gWqLLLS8Agy2e49bMGm37SUs67iMcTjSi4zpgANF0gZB88MAHgX/vNBLmqB7Wfxo40szeZ2Y7sv0lqDeJzk5EgtC/JCSPRhDdXH0X2ASc6+7/bWanAjeY2QCi//avB56L3/N3M3uSaB7qjvsG3yO6NHQh0TX7JL4db39pnAxeBrqcR9fd283su8BCovHlnwfWxS/fBdxsZucTja4pUlEafVTqnpktAC5y95Yax7GLu6+Pzwh+Adzq7oWTl4tUnC4NiaTH5WbWSjRZzl/I4DSLkk06IxARqXM6IxARqXNKBCIidU6JQESkzikRiIjUOSUCEZE69/8B8WomZYwoj+IAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(df[:50]['sepal length'], df[:50]['sepal width'], label='0')\n",
    "plt.scatter(df[50:100]['sepal length'], df[50:100]['sepal width'], label='1')\n",
    "plt.plot(test_point[0], test_point[1], 'bo', label='test_point')\n",
    "plt.xlabel('sepal length')\n",
    "plt.ylabel('sepal width')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### scikit-learn实例"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.neighbors import KNeighborsClassifier"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "KNeighborsClassifier()"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "clf_sk = KNeighborsClassifier()\n",
    "clf_sk.fit(X_train, y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "clf_sk.score(X_test, y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### sklearn.neighbors.KNeighborsClassifier\n",
    "\n",
    "- n_neighbors: 临近点个数\n",
    "- p: 距离度量\n",
    "- algorithm: 近邻算法，可选{'auto', 'ball_tree', 'kd_tree', 'brute'}\n",
    "- weights: 确定近邻的权重"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### kd树"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**kd**树是一种对k维空间中的实例点进行存储以便对其进行快速检索的树形数据结构。\n",
    "\n",
    "**kd**树是二叉树，表示对$k$维空间的一个划分（partition）。构造**kd**树相当于不断地用垂直于坐标轴的超平面将$k$维空间切分，构成一系列的k维超矩形区域。kd树的每个结点对应于一个$k$维超矩形区域。\n",
    "\n",
    "构造**kd**树的方法如下：\n",
    "\n",
    "构造根结点，使根结点对应于$k$维空间中包含所有实例点的超矩形区域；通过下面的递归方法，不断地对$k$维空间进行切分，生成子结点。在超矩形区域（结点）上选择一个坐标轴和在此坐标轴上的一个切分点，确定一个超平面，这个超平面通过选定的切分点并垂直于选定的坐标轴，将当前超矩形区域切分为左右两个子区域\n",
    "（子结点）；这时，实例被分到两个子区域。这个过程直到子区域内没有实例时终止（终止时的结点为叶结点）。在此过程中，将实例保存在相应的结点上。\n",
    "\n",
    "通常，依次选择坐标轴对空间切分，选择训练实例点在选定坐标轴上的中位数\n",
    "（median）为切分点，这样得到的**kd**树是平衡的。注意，平衡的**kd**树搜索时的效率未必是最优的。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 构造平衡kd树算法\n",
    "输入：$k$维空间数据集$T＝\\{x_1，x_2,…,x_N\\}$，\n",
    "\n",
    "其中$x_{i}=\\left(x_{i}^{(1)}, x_{i}^{(2)}, \\cdots, x_{i}^{(k)}\\right)^{\\mathrm{T}}$ ，$i＝1,2,…,N$；\n",
    "\n",
    "输出：**kd**树。\n",
    "\n",
    "（1）开始：构造根结点，根结点对应于包含$T$的$k$维空间的超矩形区域。\n",
    "\n",
    "选择$x^{(1)}$为坐标轴，以T中所有实例的$x^{(1)}$坐标的中位数为切分点，将根结点对应的超矩形区域切分为两个子区域。切分由通过切分点并与坐标轴$x^{(1)}$垂直的超平面实现。\n",
    "\n",
    "由根结点生成深度为1的左、右子结点：左子结点对应坐标$x^{(1)}$小于切分点的子区域， 右子结点对应于坐标$x^{(1)}$大于切分点的子区域。\n",
    "\n",
    "将落在切分超平面上的实例点保存在根结点。\n",
    "\n",
    "（2）重复：对深度为$j$的结点，选择$x^{(1)}$为切分的坐标轴，$l＝j(modk)+1$，以该结点的区域中所有实例的$x^{(1)}$坐标的中位数为切分点，将该结点对应的超矩形区域切分为两个子区域。切分由通过切分点并与坐标轴$x^{(1)}$垂直的超平面实现。\n",
    "\n",
    "由该结点生成深度为$j+1$的左、右子结点：左子结点对应坐标$x^{(1)}$小于切分点的子区域，右子结点对应坐标$x^{(1)}$大于切分点的子区域。\n",
    "\n",
    "将落在切分超平面上的实例点保存在该结点。\n",
    "\n",
    "（3）直到两个子区域没有实例存在时停止。从而形成**kd**树的区域划分。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# kd-tree每个结点中主要包含的数据结构如下\n",
    "class KdNode(object):\n",
    "    def __init__(self, dom_elt, split, left, right):\n",
    "        self.dom_elt = dom_elt  # k维向量节点(k维空间中的一个样本点)\n",
    "        self.split = split  # 整数（进行分割维度的序号）\n",
    "        self.left = left  # 该结点分割超平面左子空间构成的kd-tree\n",
    "        self.right = right  # 该结点分割超平面右子空间构成的kd-tree\n",
    "\n",
    "\n",
    "class KdTree(object):\n",
    "    def __init__(self, data):\n",
    "        k = len(data[0])  # 数据维度\n",
    "\n",
    "        def CreateNode(split, data_set):  # 按第split维划分数据集exset创建KdNode\n",
    "            if not data_set:  # 数据集为空\n",
    "                return None\n",
    "            # key参数的值为一个函数，此函数只有一个参数且返回一个值用来进行比较\n",
    "            # operator模块提供的itemgetter函数用于获取对象的哪些维的数据，参数为需要获取的数据在对象中的序号\n",
    "            #data_set.sort(key=itemgetter(split)) # 按要进行分割的那一维数据排序\n",
    "            data_set.sort(key=lambda x: x[split])\n",
    "            split_pos = len(data_set) // 2  # //为Python中的整数除法\n",
    "            median = data_set[split_pos]  # 中位数分割点\n",
    "            split_next = (split + 1) % k  # cycle coordinates\n",
    "\n",
    "            # 递归的创建kd树\n",
    "            return KdNode(\n",
    "                median,\n",
    "                split,\n",
    "                CreateNode(split_next, data_set[:split_pos]),  # 创建左子树\n",
    "                CreateNode(split_next, data_set[split_pos + 1:]))  # 创建右子树\n",
    "\n",
    "        self.root = CreateNode(0, data)  # 从第0维分量开始构建kd树,返回根节点\n",
    "\n",
    "\n",
    "# KDTree的前序遍历\n",
    "def preorder(root):\n",
    "    print(root.dom_elt)\n",
    "    if root.left:  # 节点不为空\n",
    "        preorder(root.left)\n",
    "    if root.right:\n",
    "        preorder(root.right)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 对构建好的kd树进行搜索，寻找与目标点最近的样本点：\n",
    "from math import sqrt\n",
    "from collections import namedtuple\n",
    "\n",
    "# 定义一个namedtuple,分别存放最近坐标点、最近距离和访问过的节点数\n",
    "result = namedtuple(\"Result_tuple\",\n",
    "                    \"nearest_point  nearest_dist  nodes_visited\")\n",
    "\n",
    "\n",
    "def find_nearest(tree, point):\n",
    "    k = len(point)  # 数据维度\n",
    "\n",
    "    def travel(kd_node, target, max_dist):\n",
    "        if kd_node is None:\n",
    "            return result([0] * k, float(\"inf\"),\n",
    "                          0)  # python中用float(\"inf\")和float(\"-inf\")表示正负无穷\n",
    "\n",
    "        nodes_visited = 1\n",
    "\n",
    "        s = kd_node.split  # 进行分割的维度\n",
    "        pivot = kd_node.dom_elt  # 进行分割的“轴”\n",
    "\n",
    "        if target[s] <= pivot[s]:  # 如果目标点第s维小于分割轴的对应值(目标离左子树更近)\n",
    "            nearer_node = kd_node.left  # 下一个访问节点为左子树根节点\n",
    "            further_node = kd_node.right  # 同时记录下右子树\n",
    "        else:  # 目标离右子树更近\n",
    "            nearer_node = kd_node.right  # 下一个访问节点为右子树根节点\n",
    "            further_node = kd_node.left\n",
    "\n",
    "        temp1 = travel(nearer_node, target, max_dist)  # 进行遍历找到包含目标点的区域\n",
    "\n",
    "        nearest = temp1.nearest_point  # 以此叶结点作为“当前最近点”\n",
    "        dist = temp1.nearest_dist  # 更新最近距离\n",
    "\n",
    "        nodes_visited += temp1.nodes_visited\n",
    "\n",
    "        if dist < max_dist:\n",
    "            max_dist = dist  # 最近点将在以目标点为球心，max_dist为半径的超球体内\n",
    "\n",
    "        temp_dist = abs(pivot[s] - target[s])  # 第s维上目标点与分割超平面的距离\n",
    "        if max_dist < temp_dist:  # 判断超球体是否与超平面相交\n",
    "            return result(nearest, dist, nodes_visited)  # 不相交则可以直接返回，不用继续判断\n",
    "\n",
    "        #----------------------------------------------------------------------\n",
    "        # 计算目标点与分割点的欧氏距离\n",
    "        temp_dist = sqrt(sum((p1 - p2)**2 for p1, p2 in zip(pivot, target)))\n",
    "\n",
    "        if temp_dist < dist:  # 如果“更近”\n",
    "            nearest = pivot  # 更新最近点\n",
    "            dist = temp_dist  # 更新最近距离\n",
    "            max_dist = dist  # 更新超球体半径\n",
    "\n",
    "        # 检查另一个子结点对应的区域是否有更近的点\n",
    "        temp2 = travel(further_node, target, max_dist)\n",
    "\n",
    "        nodes_visited += temp2.nodes_visited\n",
    "        if temp2.nearest_dist < dist:  # 如果另一个子结点内存在更近距离\n",
    "            nearest = temp2.nearest_point  # 更新最近点\n",
    "            dist = temp2.nearest_dist  # 更新最近距离\n",
    "\n",
    "        return result(nearest, dist, nodes_visited)\n",
    "\n",
    "    return travel(tree.root, point, float(\"inf\"))  # 从根节点开始递归"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 例3.2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[7, 2]\n",
      "[5, 4]\n",
      "[2, 3]\n",
      "[4, 7]\n",
      "[9, 6]\n",
      "[8, 1]\n"
     ]
    }
   ],
   "source": [
    "data = [[2,3],[5,4],[9,6],[4,7],[8,1],[7,2]]\n",
    "kd = KdTree(data)\n",
    "preorder(kd.root)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "from time import clock\n",
    "from random import random\n",
    "\n",
    "# 产生一个k维随机向量，每维分量值在0~1之间\n",
    "def random_point(k):\n",
    "    return [random() for _ in range(k)]\n",
    " \n",
    "# 产生n个k维随机向量 \n",
    "def random_points(k, n):\n",
    "    return [random_point(k) for _ in range(n)]     "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Result_tuple(nearest_point=[2, 3], nearest_dist=1.8027756377319946, nodes_visited=4)\n"
     ]
    }
   ],
   "source": [
    "ret = find_nearest(kd, [3,4.5])\n",
    "print (ret)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:2: DeprecationWarning: time.clock has been deprecated in Python 3.3 and will be removed from Python 3.8: use time.perf_counter or time.process_time instead\n",
      "  \n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "time:  5.170202400000001 s\n",
      "Result_tuple(nearest_point=[0.09916902877403755, 0.5005978535517558, 0.7997848590100571], nearest_dist=0.0010460533893058112, nodes_visited=38)\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\ProgramData\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:5: DeprecationWarning: time.clock has been deprecated in Python 3.3 and will be removed from Python 3.8: use time.perf_counter or time.process_time instead\n",
      "  \"\"\"\n"
     ]
    }
   ],
   "source": [
    "N = 400000\n",
    "t0 = clock()\n",
    "kd2 = KdTree(random_points(3, N))            # 构建包含四十万个3维空间样本点的kd树\n",
    "ret2 = find_nearest(kd2, [0.1,0.5,0.8])      # 四十万个样本点中寻找离目标最近的点\n",
    "t1 = clock()\n",
    "print (\"time: \",t1-t0, \"s\")\n",
    "print (ret2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 第3章 k近邻法-习题\n",
    "\n",
    "### 习题3.1\n",
    "&emsp;&emsp;参照图3.1，在二维空间中给出实例点，画出$k$为1和2时的$k$近邻法构成的空间划分，并对其进行比较，体会$k$值选择与模型复杂度及预测准确率的关系。\n",
    "\n",
    "**解答：**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "from sklearn.neighbors import KNeighborsClassifier\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.colors import ListedColormap\n",
    "\n",
    "data = np.array([[5, 12, 1], [6, 21, 0], [14, 5, 0], [16, 10, 0], [13, 19, 0],\n",
    "                 [13, 32, 1], [17, 27, 1], [18, 24, 1], [20, 20,\n",
    "                                                         0], [23, 14, 1],\n",
    "                 [23, 25, 1], [23, 31, 1], [26, 8, 0], [30, 17, 1],\n",
    "                 [30, 26, 1], [34, 8, 0], [34, 19, 1], [37, 28, 1]])\n",
    "X_train = data[:, 0:2]\n",
    "y_train = data[:, 2]\n",
    "\n",
    "models = (KNeighborsClassifier(n_neighbors=1, n_jobs=-1),\n",
    "          KNeighborsClassifier(n_neighbors=2, n_jobs=-1))\n",
    "models = (clf.fit(X_train, y_train) for clf in models)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2cAAAE/CAYAAADCCbvWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzde3SbV3rf++8GQRIkQRLUlZIoCbJ1oWTZA9myTc/IM8jMJOPcmqRx0s6k6aBN4/S0OeewzVmrbv4Jpj1t1HNyQdZJmmTSZGHSZJJmOUnTpMmkk8trR7bpsWzDtmxTliy9omgJkigJpCAKJEHs8wdAiDdJJEXyBcDfZy2sAV7cHtAabDzv3vt5jLUWERERERER8ZbP6wBEREREREREyZmIiIiIiEhFUHImIiIiIiJSAZSciYiIiIiIVAAlZyIiIiIiIhVAyZmIiIiIiEgFUHImFccYkzXGPLDAx1pjzO473Bczxhxb3uiWjzHmR4wx/+su90eNMYOLeD3HGPPPlic6ERGRmTQ+l+/X+CwrRsmZ3JMxxjXGfH7a7X9ojLlujPnMPI+Nlr6Qf2XW8WPGmNhC3s9aG7TWnrnvwCuctfZ3rbXfMXX7bgPZajHG/LAx5hVjzKgxxvEyFhERuTuNzyujQsfnnzPGnDLG3DDG9Btj/rGX8cjKUXImi2KM+TLwK8B3W2tfvMPDbgL/2BgTXq24Vpoxxu91DKvkGpAAjnodiIiILJzG55p3E/heoB34MvBLxphPehuSrAQlZ7JgxpjngJ8HvmCtfeUuD80ASeBn7vJa/9QY80HpDN9fGmN2TruvfIbKGLPeGPOnxpgRY8zrxpj/e56lEJ8vnU26boz5FWOMmflW5v8zxgyXzjR9btodW40x/8MYc80Yc9oY8+PT7osbY14wxvyOMWYEiBljnjDGHC/FcskY8wt3+GwvGmN+sHT9SOnzfFfp9ueNManS9fKyDmPMS6Wnv11aNvIPpr3eTxljLhtjLhpj/skd/+ozY9hijHnHGPN/LeTxU6y1f2Wt/QPgwmKeJyIi3tH4vCbG55+x1vZbawvW2teAvwOeWsxrSHVQciYL9b8B/x74nLX2+AIe/x+AHzTG7Jt9hzHm+4GfBv4+sJHiF8zv3eF1foXi2aJOimeKvjzPY74HeBz4BPDDwBem3fckcAbYQHEw+iNjzLrSfb8HDAJbgWeB/zh9cAC+D3gBCAG/C/wS8EvW2jbgQeAP7hDzi0C0dP3Tpff/zLTbc85oWms/Xbr6idKykf9Wut1J8SzZNuDHgF8xxnTc4X0BKJ0RfRH4ZWvtz5WO/WdjTOYOl3fu9noiIlLRND6vsfHZGNNE8e/63t3eT6qTkjNZqG8H+oB3F/Jga20a+DXg381z908AP2ut/cBamwf+IxCZfnYOwBhTB/wg8DPW2lFr7fvA1+Z5vaPW2oy1dgD4WyAy7b7LQMJaO1H6Qj0JfLcxZjtwBPg31tqctTYF/BfgR6c991Vr7X8vnaW6BUwAu40xG6y1WWtt3x0+/ovM/LL/2Wm3P8M8X/53MQH8u1L8fw5kgTkD6jQHAIfi3+yrUwettf/CWhu6w+WRRcQjIiKVRePz2huffw14G/jLRcQrVULJmSzUPwf2Av9l1rKEu/lPwBeMMZ+YdXwnxbXSGWNMhuI+J0Px7NN0GwE/cH7asfPMlZ52fRQITrv9sbXWTrt9juKZuK3ANWvtjVn3TY9h9nv9GMW/QX9pCcf3zBMLwKvAXmPMZooD0W8D240xG4AngJfu8Lz5XC0NkFNmf77ZfgT4mOIZRRERqX0an9fQ+GyM+X+Bg8APz/r7SY1QciYLdRn4HPA08J8X8gRr7VWKxSX+/ay7zgM/MevsUNM86+SvAHmga9qx7YuMe9uswWoHxf1UF4B1xpjWWfd9PP0jzPo8p6y1XwQ2URzYXjDGtMx+Q2vtKPAG8H8CJ6y148ArwL8GPrLWDi3yMyxGHBgCvl46swmAMebXSmvl57toWYSISPXS+LxGxmdjzFeA7wS+w1o7soKxioeUnMmCWWsvAJ8FnjHG/OICn/YLwCeB/dOO/Rrwb40xDwEYY9qNMT80z/tNAn8ExI0xzcaYbmCxpWM3Af+HMaa+9B77gT+31p6n+IX8s8aYgDHmEYpn3n73Ti9kjPlHxpiN1toCxU3VAJN3ePiLwE9ye4mEM+v2fC4BC+ofcxcTwA8BLcB/Ncb4AKy1/7y0Vn6+y0NTTzbG1BljAhTPiPpKf5v6+4xJRERWkMbnNTE+/1vgS8C3l5JrqVFKzmRRSl+anwWeNcb87AIePwL8P8C6acf+mOKZrd83xUpLJyieCZrPT1LccJsG/ivFTcJjiwj5NWAPxbNV/wF4dtqX2heBMMWzdH9McR34N+/yWs8A7xljshQ3H/9Da23uDo99EWjl9hKJ2bfnEwe+VlpO8sP3+Fx3VDoT+PcpDny/NTUALNCPAreAX6V4FvYW8BtLjUVERFaHxueaH5//I8UZxFPTZtZ+eqmxSOUyWq4q1cQY85+ATmvtfFWhRERExAMan0WWh2bOpKIZY7qNMY+YoicoLm34Y6/jEhERWcs0PousjLXSVV2qVyvFpRJbKW56/nngTzyNSERERDQ+i6wALWsUERERERGpAFrWKCIiIiIiUgGUnImIiIiIiFSAVd1z1tzebEOdodV8SxGRNevihxeHrLUbvY5DKt+G5mYbDq3S+Hz1Kql141BXx6aWTavzniIiFeRu4/OqJmehzhDP/fpzq/mWIiJr1le+7SvnvI5BqkM4FOL4c6s4PieThJ8dJLNhlN6e3tV7XxGRCnC38VnLGkVERGR1xWK4L3RB7k59gkVE1iYlZzXIWsvwpWEG3h3g4w8+ZvzWuNchiYiIzMtxHa9DWFW3btzi/HvnGXh3gOy1rNfhiEiFUZ+zGpMfz/PGN9/g4tBF7HqLyRt8r/iI9ETY8dAOr8MTEREpCocJZQdJuX0ARMNRb+NZYdZaTr1+ig/e/QC70YIBjsOu8C4e/vTD+Op0vlxElJzVnPdffp8LExcIfTqEMQaA/GieN7/1Jq0drXRs7fA4QhERESAaxXUg2t1HitpP0C5+eJETH56g/el26hrqALCTlo+Of0RLqoXdj+32OEIRqQQ6TVNDxm+N455xaX+ovZyYAfib/dSF6zhz4oyH0YmIiMwSjeJ8o5MQAa8jWXGn3j1F096mcmIGYOoMrQ+1curEKWzBehidiFQKJWc1JJfNQQB89XP/szaua2T4+rAHUYmIiNxDLoebcb2OYkWNZEZoXNc453h9sJ7x/DgTYxMeRCUilUbJWQ1pbGnE5iyFfGHOfePD4wTbgh5EJSIichexGOEMZIYGSaaSXkezYlpaWxgfnlugKz+ax+/z42/QThMRUXJWUxqbG+na3sXIhyNYe3t5xOT4JBNnJ3jgoQc8jE5ERGR+Tufz9J4Iksmka7Z6456Dexj9cBQ7eXt8tgXLSP8ID3Y/qIIgIgKoIEjNefjph7n5P29yte8qvg0+7ISFS3Dw4EHWb1/vdXgiIiLziqdCJCMZr8NYMV37u7h2+Rpnj52FzWB8BnvJsm3dNvYc3uN1eCJSIZSc1ZiGpgaO/MARrp6/ytWLV6lvqKfzqU5aOlq8Dk1EROTusllS6VRNVm00PsMj0UcIXwlz2b2MtZYND22gY2vHjCJeIrK2aQ69BvnqfGwMb6T7qW4efOxBJWYVqDBZoDA5d2+giMiaFYsROx2EkRESfQmvo1kRxhjaN7Wz54k97H1yL+u2rVNiVmFswWp8Fk9p5kxkFY1cGeHk8ZNcOH8BgK07ttJ9uJvWDa0eRyYi4r14pJd4MkkoVtx7VoszaFKZxm+N8+HrH+J+6JKfzNOxvoPuR7vZ/MBmr0OTNUYzZyKrZOTKCC/92Uukm9K0fraV1s+2kg6keelPX+LG1RtehyciUjHWQt8zqRz58Twv/4+XOX3jNE1PN9H+He2M7hjllRdfYbB/0OvwZI1RciaySk6+cRK7y9IabsXn9+Hz+2gNt1IIF/jwjQ+9Dk9EpHJks6Tcvpqt3CiV5eOTH5OpzxB6KERdYx3GGJo2NdHyWAvvvvauljnKqlJyJrIKCpMFLgxcoKVr7v6/lu0tDJ4dnNH+QERkzYrFcF/oIpTNK0GTVTF4dpBA19zZ2oa2Bsbrxxm5MuJBVLJWKTkTWQXa8C0isghTCZqWN4rXdN5UVpmSM5FVYHyGbTu2cfP8zTn3ZQeydO3qUgInIjJbLud1BLIGbN+1ndz5uf/WxofHaZxspG1jmwdRyVql5Exklex7fB8+18fImREK+QKFfIGRMyPUnatj3+F9XocnIlJZwmHCGUidPkYylfQ6GqlhW/dtpWOyg8y7GfK38tiCZTQ9ys03b/LIk4/gq9PPZVk9+tcmskpa17fyme/7DNvy28j+TZbs32TpmuziM9/3GYLrgl6HJyJSWaJRnP4eIkN+MkODStBkxfgb/Hzyez/JvnX7GHt1jOFvDtN2sY2nP/c0W/dt9To8WWPU50xkFQXXBXnsOx7jUfsooL1oIiJ3FY3iECWeSpA4WEzQYpGY11FJDWpoamD/p/bT/clusMXtCCJe0MyZiAeMMUrMREQWKB7ppfdEUDNosuKMMUrMxFP3TM6MMQFjzLeMMW8bY94zxnyldHyXMeY1Y8wpY8x/M8Y0rHy4IiIiAmtvfI5nIkSyWgIuIrVtITNnY8BnrbWfACLAM8aYHuA/Ab9ord0DXAd+bOXCFBERkVnW3vicy5HJZbyOQkRkxdwzObNF2dLN+tLFAp8FXigd/xrw/SsSoYiIiMyx5sbnaBTnG50wMkKiL+F1NCIiK2JBe86MMXXGmBRwGfgm8BGQsdbmSw8ZBLbd4bnPGWOOG2OOX75ymWQqqfXiIiIiy2C5xucro6OrE/D9isXIfH0HoUxOvyVEpCYtKDmz1k5aayNAF/AEsH++h93huV+11h621h5uqoewm9GGXhERkWWwXOPzxubmlQxzeYXDhHMBr6MQEVkRi6rWaK3NAA7QA4SMMVOl+LuAC/d6/r7RJpxUZEbFJc2kiYiI3J/7HZ+rTi5HJpPGcR2vIxERWVYLqda40RgTKl1vAj4PfAD8LfBs6WFfBv7knu/W2grRaLkkbrg/Tbg/TebygNaPi4iILMKyjs/VpNScOpTNk3L7lKCJSE1ZyMzZFuBvjTHvAK8D37TW/hnwb4B/bYw5DawHfnMxbxyP9OJ0Po/T+Ty977dpg6+IiMjirMj4XBWiUdzjR4hkArgZ1+toRESWjf9eD7DWvgMcmuf4GYrr2+9bPNJLPJkk9KXiDFpvT+9yvKyIiEjNWo3xWUREVtc9k7NVE4uRSSYJPztIwjkKgeJm30hnhGg46m1sIiIiUnlyufIe9lgk5nU0IiL3bVEFQVZcLIb7Qhe9qQC9fRAaymo9uYiIiMwVjeJ0Pk9kyK8q0CJSMypn5mxKLEa8dDXuOIQPHyPl9gFoBk1ERERmcPp7iEZSuBu8jkRE5P5V1szZbFMbftOQOn1MM2giIiIyVy7ndQQiIsuispMzKJfMjQz5SZ0+RsI5SsI5quULIiIiAtEo0XRAbXlEpCZUfnIG5XXlmd/pIpPspDdV/BJWgiYiIiLxSG+5LY9+G4hINauO5GxKLFbck1b6EtZZMhEREYFigha51uB1GCIi96W6krNppp8lS/QlcFxHe9JERETWuEwu43UIIiJLVrXJGRQTtMzXdxAayuL29xX3pGkmTUREZE1yvtFJaChb7JcqIlKFKq+U/mLFYriOU7zuuuUm1qFQJ+FQGFAJfhGpHSNXRjjff54bIzdo72hnx/4dtHS0eB2WSGWIxXCTyfJvgUi4R78BRGRV5LI5BvsHuXr5KoGmANv3bqdjawfGmEW9TvUnZwDRaPmqm0wSfSYNmQxuoI9MoPQQfTmLSJUbeG+AN7/1Jr4uH/71fi5lLnHqj0/xZPRJNj+w2evwRCpD6aSt+qSKyGoZvjTMsT8/Rn5Tnvr19eRH85z9q7N07+1mX8++RSVotZGcTReL4UzNpAHR7j5S+WPF6/pyFpEqdWvkFm+99hbBniD+5tJXdyeMbxnnded1nul6Bn9D7X2liyxJNIrrlH4DBFIa/0VkxdiC5fhfH8fsN7R3tpePF7YX6H+ln807N9OxtWPBr1fVe87uKBotX6b3SFN5XRGpVhc/uojdZG8nZiUN7Q3k2/Ncca94FJlIhSr1PxMRWUnDl4e5YW/Q3Nk847iv3odvu4+BkwOLer3aTM6mm5agZYYGSaaS5YuISLUYuzWGr+kOX9kBmBibWN2ARKpFNqsxX0RWzERuAhOYf9miv9lP7lZuUa9X+8kZlJtY954IEu5PE+5Pq4m1iFSVjk0dFK4V5hy31sI1aF3f6kFUIpUtHumdcXJWRGS5BdcHsRlLIT93jB6/Ms76zesX9XprIzkriUd6cTqfLyZqamItIlVkY3gjrROt3Dhzo5iQUVznPnJyhA0tGwhtCXkcoUhlmjo5q/5nIrISmlqbCIfDDJ8YLido1lpG06M0XG2ga1/Xol5vTSVn081uYi0iUsnq/HU89d1PEcqEGH5pmOE3hhl+cZjN+c08/szjiy7VKyIiIsvj4U8/zK7WXYy8OMLIGyOMvDxC49lGPvVdnyIQXNze1zVd2ise6SWeTBL60swZtN6eXg+jEhGZX3N7M0d+4AjZq1ly2RzN7c3qcSayUKWTsRrjRWS51dXXEflchH039pG9lqW+sZ72ze1LOnG6ZmfOymIxMl/fQW8f9PZBaChLwjmK4zpeR1bzCpMFMukM1z6+Rn4873U4IlXBGEPrhlY2hjcqMRNZIK2WWRxrLTeGbnB18CpjN8e8DkekajS1NrFx50ZCnaElr2hZ0zNnZbEY8dLVeDJJ+NlBNa5cYemP0qReTjFWNwZ14L/l58ChA4Q/EdbyLBERWXbxSC+kEiR6vI6ksmWvZXnjb97g+s3rmCYDNyD8QJiDRw5SV1/ndXgiNU/J2WyxGK7jED58TAnaCrl+4Tp9L/XRdKiJ9lCxWV9+NM9bb7yFv8HP9gPbPY5QRERk7Rm/Nc7Lf/YyE7smaO8qLskq5AucefcMhRcLHPr8Ia9DFKl5WtY4n2gU9/gRImlInT5GwjlKwjmqMrzL5FTqFP4H/TSGGsvH/M1+Wh5u4YM3P8AWrIfRiYhIrYqnQtq+cBcXPrzArbZbBLcHy6tYfH4foUdCnDt3jtHhUY8jFKl9Ss7upNS8OvM7XWSSnUTSlPukOK6jL/X7cCV9habNTXOON4YaGR0bZTw37kFUIiJS82Ix3Be6CGXzpNIpr6OpOEOXhqjfWD/nuKkzmA7Djas3PIhKZG3Rssa7iUbLVx0gnkqQ3D2Im8mQIUcqnVLVpyVoDDQymZukrnHm2vXCRAFfwUedX2vaRURkhcRixFIJEhu8DqTyTI3P8xoDf4N+NoqsNM2cLUI80ov7QhduMkRvKqCqT0u0a98ubp6+WW6kO+XG2Rt0hbv05S8iIisvl9MqmFm69nRRGCxQmCjMOJ4byhGYCNCxpcOjyETWDiVnixWLFas7qizvku18eCeb6zaT+VaG0Yuj3Lp0i0wqQ3AoyIGnDngdnoiI1Lh4JlLeV6795LeFOkPs27OP4VeHuXHuBrmhHMP9w0y8O8Hjn3scX51+NoqsNE1R3IfZTaxDgRAAsUjM28AqnL/Bz5Pf8yTpj9KcP32eQqHA1ge2sm3vNuoDc9e6i4iILKtoFMeBaHcfKX/a62gqhjGG7qe62bR9E+dOnuPWx7cIbwqz46kdNLc3ex2eyJqg5Ox+xWJkkkmiz6QhkMEN5Eg4R4mEe1SC/y7q/HVs27eNbfu2eR2KiIisRdEoTtIlFFNyNp0xhvXb17N++3qvQxFZk5ScLYdYDMdxyjfVI01kaQqTBT7u/xj3Q5fxsXE2b93Mrod30dLR4nVoIlKr8nmSqaRWvYjchbWWoYEhzpw4Q/ZGlrZQGw8cfID1XUril5sWDy+XaLR8md4jTZuNRRamMFng9b94nddPvE52W5bJA5N8lPuIv/3vf8v1i9e9Dk9EalEsRu+JIJnLA9o/LnIXp4+f5tiLxxhqH2LywCSXmi/x0l+9xNm3z3odWs1RcrYSSj3SIkP+8mbjqYuIzO/iqYt8fPNjOp7ooGlTEw1tDbTvbce330fqpdSc6p4iIstBBb5E7u5m5ibvv/s+bT1ttHS10NDWQHBHkNYnW3n3+LuM3RzzOsSaouRspUxL0ML9acL96XITaxGZ69ypcwR2BjDGzDjetLmJ4dwwN6/f9CgyEal1StBE7uzS2UvYzZa6hpl9aP1NfgobClw5d8WjyGqTkrOVFI3idD5fvkwtnVCCJjLXxPjEnC9+KG5ONw2GyYk7NEYVEVkG8Ugvma/vUP8zkVnyE3lMvZn3PltvNT4vMyVnq2jqzJzWtovM1dnVSe5ibs7x/Ggef86voiAisipCBLwOQaSirOtch71i52wvsNZihgztm9s9iqw2KTlbZVo6ITK/nQ/tpOFKA9mBLLZQHAAmbk4w8tYI3Z/oxt+g4rIisgqyWVJuH4m+BIm+hGbRZM1b37WeDc0bGH5vmEK+AEBhokDmnQxb1m9RcrbMlJx5oLx0QgmaSFkgGODI9x4hdDXE8IvDDB8bZuL1CSIPRXjg0ANehycia0EshvtCF72pAL19EBoqJmpK0GQtMz7Dk9/5JDubdpJ9Mcvwy8NkX8ryYPuDPPbtj83ZKy73R6eivRKL0ZtKkOjxOhCRytG6vpUjP3CEWyO3yE/kaW5vps4/dx+aiMiKicWIl67GHYdodx+p/O3+pVNCoU71RpM1oz5Qz6HPHeKhWw8xdnOMQDBAfaDe67BqkpIzr+Xm7rGR2jZ2c4xrH18DYN22dTS2NHocUeVpamvyOgQRkWJhLwdw3ZmHn0mTyg+qeXWNmZyY5OrgVfLjedo2thFcF/Q6pIrT0NRAQ1OD12HUtHsmZ8aY7cBvA51AAfiqtfaXjDFx4MeBqfqZP22t/fOVCrQWxVMhkrsHSThHiYR7iIajXockK8hay6nXT9H/bj92ncVi8f2dj/2P7Gf34d1aFiAii6LxeZVEo3MOOUA8lSBxYEAJWo24dOYSx188zkRwAhqAV6BrSxeRz0a051lW1UL2nOWBn7LW7gd6gH9pjDlQuu8XrbWR0kVf/IsVi+EeP0Iom5+zXEJqz8cffMyJD0/Q8nQL7YfaCR0K0XKkhXdPvsuFkxe8Dk9Eqo/GZw+pAnPtyF7L8tqLr1H/aD2hx0OEPhGi/TPtDOYGOXHshNfhyRpzz+TMWnvRWvtm6foN4ANg20oHtmZEo7gvdHkdhawway0n3z5Jy4GWGb286hrraN7fzMm3T3oYnYhUI43P3lMF5tpw7r1z2C5LQ/vt5XrGZ2g72Ma5s+cYuznmYXSy1iyqWqMxJgwcAl4rHfpJY8w7xpjfMsZ0LHNsa0s+j+M65YvUFluw3Bi5QUPH3HXajesaGbk+Ui4fLyKyWBqfvaMKzNXv2tA1GtfN3f/t8/sgCKPDox5EJWvVgpMzY0wQ+EOg11o7Avwq8CAQAS4CP3+H5z1njDlujDl+ZVT/uOcVDhMZ8uP29+H295E6fYxkKul1VLKMjM/Q0NhAfjQ/5778zTyBpgDGpz1nIrJ4Gp8rQCxG5us7CA1lSThHvY5GFqkl2MJEdmLOcVuw2FFLQ7MKYMjqWVByZoypp/jF/7vW2j8CsNZestZOWmsLwG8AT8z3XGvtV621h621hzc2Ny9X3LUlGsXpfB63rwe3r4feE0Eyl4ubjKcuUt2MMew+sJsbJ29g7e0ZMluwZE9mefDAgx5GJyLVSuNzBSn1SAtl8ySco1oFU0V2du8k7+YpTBRmHM+6WTZ2bKQl1OJRZLIWLaRaowF+E/jAWvsL045vsdZeLN38AUA7Ju9XqSJUnCikEjjX0hAIkApmSfQl6O3p9TQ8uT8PHnqQa5eucfHVi9R1FvedFS4W2BraqibLIrJoGp8rUCyG6ziED9/ui6ZKzJVvXdc6Htr/EO8fex+z1eAL+JgcmiSYC3Loew55HZ6sMQupDfop4EeBd40xqdKxnwa+aIyJABZwgZ9YkQjXqHikl2JzFYiHUiQOFNeyhwIhAJXtrUJ19XU8+d1PcnXwKhfdi/iMj81Pb2Z913otaRSRpdD4XImiUVyHcvNqUIJW6Ywx7H1iL527Orlw+gLjY+OsO7COzgc771pGX20UZCWY6UusVtrhrVvt8eeeW7X3qynJJNFn0gC4wTyZoJ/e6PMeByUilewr3/aVN6y1h72OQyqfxucV4DjFBG1DntCGLv2IrzGJvgRks+D3q1etLNrdxmd11asWsRjO1PVkkvCzxebVoVAnAOFQuCa/GAqTBQZODHD6vdOMZkdpX9fO3k/spXN3p5o2i4hI5YpGcZzSDBqDNTnLcvnsZT5Mfci1K9cINAfYfWA3Ox/ZSZ2/7t5PrjKO6+BmXAAymTShbB73hS7Czw5qCassq0WV0pcKUWpeHUlDuD8Ng8UvhlrbfGyt5a2/eos3+99k8sAkwc8GGd0xyquvvsqZN894HZ6IiMjdlQp+RYb8ZIYGa6rA17l3z/HySy8zsmWE4GeD8DCkzqQ4/pfHKUwW7v0CVcRxnWICNjhIuD9NJE2xR23p91gom6/J32HiDc2cVatoFIdo8frU0on8MVLp4raDUCBU9Wforl+4zvnL5+n4VEd5T1bTpiYa2hs48Xcn6NrfRWPz3L4kIiIilcTpfJ54KkHiwEBNzKBNjE3wzrfeofWpVvzNxZ+SDe0NdDzawYW+C1w9f5WN4Y0eR3l/kqkkmVymeCObJTLkx+k/Ui7eRqz0wKk9hpFUeWZN5H5o5qwWRKM4/cUS/L19EHFzNXGGLu2mMVvMnGIZdY11sB6unr/qUWQiIiKLE4/00vt+G5nLA1XfrPr6hesU2grlxGyK8Rn8W/1cOHvBo8iWRzKVJDM0SMTN0dsHvSeCOP09txMzkRWkmbNaEY0WS/CX1MIZOmvtHasYWrN6hWxERESWQzzSC6lEuQLzlGprlWOtvePpfR+L/mcAACAASURBVOMzc/qFVYN5Z8o6n4fOBb5A7vaJ8Wr93SWVQTNnNaoWztBt2r6JyfQksyuKFiYKmKuGjq0dHkUmIiKyNPFIL5mv7yjOyPRBaChLwjnqdViL0rGlA5MxTI5NzjhurSV/Mc+WnVs8imxpEn0JMpcHyv9Nek8Ei4nZQtXw3kJZfZo5q2Gzz9BV25m5Dds30BnsJP12mtZ9rfib/IyPjJN9L8u+7n00tTZ5HaKIiMjixWLES1fj0yowT6n00uwNTQ10P9zNieMnaHmohcZQI5O5SUZOjbDBv4FNuzbNeHyiLwG53D1fdzU/94yYxsfpfb+t+LvpPtTa3kLxhpKzGlfNCZrxGR5/5nFOHT/FmVfPkC1kaWps4tFHHmXnwzu9Dk9EROT+xWK4jgOuC1A1pdn3PL6HQEuA/rf6yeQy+I2fvfv2svfxvfjqiguzpqocTpWdv5vV/NyJvgSMjJD5+o7bB2OxZXnt8u8ulKDJ0ig5WwPikd7bZ+aqLEHzN/jZ/8n97HtyH5P5Sfz1/jvuQxMREalK0wpNuNMrMJeSlVCos+J+5Btj2PHQDrYf2E5+PE+dvw5fne922fmSUDaPe/wIxKJ3fT3XcQgfPrYiCdq8Mb2wY9kSstnikV6c9FHcTfd+rMhsSs7WilgMd9rSiUpfMjGbr85XPhMnIiJSs0rNq6dm0qLPpEnli/uYwqFw8VgFjd/GGF6++HL5dur0sWIxjW+UKmmEwwurcjhVkr6UmM64awmfd3rPsXln71YoMZuuXGBEZBGUnK0l0xK0algyISIisiZNS2YcihWYk7sHcTMZMuRwM27FzKQlU0kymTQhApDLEcn4i2Xn7zFTNq9SYhrt7sPN90EgsKTPOzcmij3KlhLTEjnf6KzKFUviPTO7Et5KOrx1qz3+3HOr9n5yB6WlA5mgv+pm0ERk4b7ybV95w1p72Os4pPJpfK4CySQA8UiGxIERaGsj0hkBVu9E6/TZKAA34xarHL7fRjwVKh5crhmpWZ83tGnHXWcOp2JzM26xR9lSZu+WW+mEeGZDUAmazHC38VkzZ2vRPEsHlKCJiIhUsFLSEwdIJUjuzuJm+sj486TSqRX/8Z/oS0A2Syg/7adjPk/v6VKVw8gyv+Hszzs+iJtOz/t5p5pGT8UWyZR6lMWWOabFisWIpRIkNngch1QVJWdr1bSlAymOVdQSiUriuI4SVxERqSjxSC9xxynecN3yfvJQqDhTFA6Fl2XsmurXlcllblc3DIdnPmgVZqSmf954KFWuQB0KhMqx9b7fRjwTWbWYFiWX0+8JWTAlZ2tZNIpDlGj6KCnU1X66cmWn/OqckRQREVmU6RUek0miz6Qhk8EN5Ehl0sWHLDEZKJfAz0E4V9y35Xxj5aobLkjp88aJEp/6vIEM5HJE0/ffo2ylxDMRnHQfKbTXXxZGyZmoaWKJ4zq4GReguJE4myd2OjjjDB2wZv8+IiJSoWIxnKmZNCiXpJ8a0xYykzbfGOi+0HV7pmwVi2nc06zPW3EzZdNNrVSiuJVEK5XkXpScCTCzaeJarCw0tV49MnT7/xJO/5HiF34qgdOZA9Kk1o2vyb+PiIhUuBm90orbFkgXZ9BSGwbvmhQ4rlMugV8+tsrVDRetkhOy2WZsJdFKJbk7JWdSVk7QSjNFtZ6ATK2lB25Xdurvuf2FXyryNGOpxLReccu9tl9ERGRZlLYtlG9O274wn7uNgbJMpm8lCar/mdyZuvrKDPFIb3HD70gxQatVyVSSzOUBwv1pwv1pek8Ei5Wd7nUmLhbDfaGLSBrC/WkYLPaMm11eWEREpFI4nc/TeyJYHvNmXxY8Bsp9i6YDXocgFU4zZzJXLEZm2gxRb/R5ryNaFjOSzanKTlOzYos5QxiL4Uxddxy1JBARkYp314IZmiVbXSNrY4WSLI1mzmR+pRmiUDZPwjla9TNDib5EMSHro3h5f5kqO0WjOP09RIb8pE4fu+OSEREREZF4pJfe99tqfoWSLJ1mzuTOYjHc0gxayq2+ErCzZ8oyX1+hMsCzNvpOvW8oENKGXxEREZmhvMe/x+tIpBIpOZO7i8Vwp5buVUmPjnJ/llIpfIB4aoX7s5Q2+sZTtxPCtd6aQEREREQWR8mZ3FsV9OgoN40uCWXzuMeP3N7cHFmdOGYslZxqTeAcLd4OBLS+XERERIinQiR319beflke2nMmCzNtb1Vm6M7leL2QTCWL/VnSkEl2kkl2zkzMPDJV+TKT7KQ3FdD6chERESmqsb39snw0cyYLN71Hh8dNFJOpJJlMmhAByOWIZEr9WSqtYWZpKWUciE9VwFSFJhEREanyvf2yMjRzJos21S8lc3lgVWfQHNfBcZ1yj7LeVAA3GcL9/c7q6M8ydZZsKFs+SzZ1ERERkTUoFsM9foRImhnbM2Tt0syZLEm50hADqzITlOhLQC5XnCnLZuk9XSqFv0p7yZbNtLNkrpsCIEMO0NkyERGRNSkaxUm6hGJpryORClCVydnI2BhvfvwxV4aHaWtu5lBXF5taWrwOa80pJ2gHVqaZ4tSMUiqdut00OhUCQitbeXGllSpg4roARJ9Jl5tYT1GiJiLVaHxyknfSadzLl2loaOChLVt4oKMDY4zXoYmIVIWqS84+unaNP3z1VQ5OTNDt93NlcpKv9ffzbY89xuGuLq/DW3PikV7iySShLy3vDFrCOUqoOKFEKJ/HfaFUCr/aZsruZNoSTKfUqsDNF5czZPz5iqyIKSJyN8O5HMlXXmHzyAj7/X5uFQp849QptjzwAN//8MP4lKCJiNxTVSVnE5OT/NG3vsU/qKtj57SZskMTE/zGG2/wwPr1rGtq8jDCNSoWIzNV7GKJJWEd18HNuADFQh/ZPO4LXRAOl94jumzhVpxSoRUcp3iz1MxaPdJEpJr82bvv8uiNGzwdCpWPPVYo8NunT/P2xo0c2rLFw+hEKpjjEH1GSxqlqKoKgpy8epUtuRw7A4EZxzvq64lYS+rCBY8ik/spCTvVoyzsZgi7GSJpiqXwY7HiDFOlF/pYLqXP6lXBFRGRpboxNsbghQv0tLbOOF7v8/F0IMBbZ896FJlIhSutnkltyBMKdXodjVSAqpo5y46Ps/4O962vq+PjW7dWNR6ZpbSXKnz42D1Lwk5POjJDg0SG/Dj909YsrpWE7A5mF1wJBYpnojWTJiKV6ObEBG0Uk7HZ1tfXkx0dXf2gRKqB65KKQiR8RPvNBaiy5GxzSwtvAdbaOZuL3clJutrbvQlMbotGcZ3S0rxSkYvZXzbJVLKckBWVepSt8YRstqkEzenMAWlS68bVI01EKlJHIMBwXR3ZyUmCdXUz7nNzOTZv2+ZRZCJVIBBQYiZlVZWchUMhzPr1HLt2jSNtbeUE7YObNznb1MR3az17ZYhGcZypvVPHynvJpmQuDxQrL0amJRmayZ/XjL/RtH19U0sfwqGwvtBFxHONfj+R3bv50w8+4Nn29vIM2tDEBC8WCvzgAw94HKGISHWoquTMGMMXH3+cP3jjDVJXrrDdGK5Yy83WVr70xBME/FX1cWpbKUGLh1I4nbc3uUbTAWBWYiYLU+qRFo9kcDrTuME8qUzxb6sETUS89vm9e/nT8XF+8cwZdhvDLWDQ7+cLTzzBzmlFQkRE5M6qLptpDwT4Z5/8JIMjIwyNjnKwsZEHOjpUorcSRaPEic48phmy+xOLEZ+6vsD9fSIiq8Hv8/EDjzzCtT17GBgept7n49l162jUiVOR+ZVWxEDQ60ikglTlN6Yxhu3t7WzXHjNZooK1XMpmsRT3MtbNs4m94i1gf5+IyGpb19SktjZyX66OjnIrn2dDc3Ptrooq9YilrU17yWWGe/6LN8ZsB36b4pxHAfiqtfaXjDHrgP8GhAEX+GFr7fWVC1Vkebx/+TL/6+238d+8iQ+41dTEZx9+mENbt3od2uLN2t+XSqcACAVCquwoUuM0PkutuZTN8qdvv83wlSu0GcNVn49De/fy+T17qvMk6l3EIxklZjKvhfxLzwM/Za3dD/QA/9IYcwB4Hvhra+0e4K9Lt0Uq2pnr1/mLV17hBycn+clQiH8RCvEjxvDia6/x3qVLXoe3NNN6o/X2QW8f6pEmsjZofJaakR0f57++/DKPXrvGv2pv58fb2/nJ5mauvP8+3+jv9zo8kVVzz+TMWnvRWvtm6foN4ANgG/B9wNdKD/sa8P0rFaTIcnnp5Em+UF/P9mmNzDsbGvh7gQAv9fdjrfUwuvsTj/SWL73vt5G5XOyRJiK1SeOz1JLjg4N037rFo62t5ToCwbo6nm1v593Tp8mOj3sc4fKKp0KEhrIknKNehyIVZlFzxMaYMHAIeA3YbK29CMUBAti03MGJLLfzly+zb569ELsCAa5dv8745KQHUS2/qQSNkRHNoImsARqfpdqdv3yZfQ0Nc44HfD62W8vHIyMeRLWCYjHcF7oIZfMknKM4ruN1RFIhFpycGWOCwB8CvdbaBf8/xBjznDHmuDHm+JXR0aXEKLJsGvx+RguFOcfHrAWfr3bWtDtOsXm1z0c4FPY6GhFZQRqfpRY0NDTMOz4D3AQaZjU3rwnTErSU26cETYAFJmfGmHqKX/y/a639o9LhS8aYLaX7twCX53uutfar1trD1trDG5ublyNmkSU7uGsXr2Szc46/duMG+3bswF8jyVm0u49UJ0R2H1EFR5EapvFZasXBri5em5hgctb2gnO5HNmmptrtlReL4R4/QiQNqdPHcFynfJG16Z6/RI0xBvhN4ANr7S9Mu+t/AF8uXf8y8CfLH57I8oo++CCn29v54+vXOZfLcT6X439mMrzZ0sLnu7u9Dm9ZhUKdSsxEapjGZ6kl+zduJBQO87VMhpOjo6THxzk2PMwfTEzwvY89Vtv9bKNRnP4eIkN+3P4+XDdFyu3TtoQ1aiHNIz4F/CjwrjEmVTr208BR4A+MMT8GDAA/tDIhiiyfloYG/tmRI7w+OMhfnj+PtZY9u3fz4zt2EJxnrbuISAXT+Cw1w2cMz0YivLNlC6+ePcutsTG2bt3KPw6H2RxcA02aS61xcN3izWfSpEjjuI5OtK4x90zOrLXHgDudrvjc8oYjsvKa6uv59K5dfHrXLq9DERFZMo3PUmt8xhDp7CTS2el1KN6IRstXnWSScCzjXSzimRptu762DedyvDYwwMClSzTU13Nw504+sXlz7RS7EPFYJp3Bfc9lZGSEtvY2wgfChDprdD+EiCyb8clJ3rxwgQ/On6dQKLB72zYe7+qiub7e69CkUjgOUGxSnSHnbSxVaHR4lIH3B7iSvkJDYwPhfWE27dqE8VXPslglZzXm4o0b/M6xY3xibIxvDwS4VSjQd/Ei73V18cXHHquZghciXjn37jneOv4Wvp0+Gnc2kslkcP/C5dEnHmXHQzu8Dk9EKlQunyf56qt0DA3xdCBAHfDu22/z6x99xD85coTQtP6bsgY5DtHuPtye4s2MP08o1KUljYuQSWd4+S9eJt+Zp3FHIyO5ES68doGwGybybZGqSdCUnNWYP3v7bb5jcpJPTKtqtLepid8ZHOStLVt4vKvLw+hkNcRTCVIH82geZ/nlsjne/tbbBD8ZxN9U/PpsXNdIvjNP6tUUm8ObaWxp9DhKEalEf3f2LFuGhvh7oRCmVNxiV1MTLw4P883+fn4oEvE4Qll1jjNzj9mGPJHdR8p3KzFbOGstbzpvYvYb2jvby8ebtzTjvuqy7dw2Nu2qjpaPSs5qyNXRUUaGhni4vX3GcZ8xfDIQ4EXXVXJW4+KpBIkDIxBsIxaJeR1Ozbl05hKFjYVyYjbF3+ynsKHA5bOX2X5wu0fRiUgle+fMGb4cDJYTsylPtrbyCwMDTDz8MPW12MtL5pdMEn52EKLFYicZIBJW+5ulujF0g5GJEdo3z/wNbHyG+h31DJwcUHImqy+XzxM0Zt5ys8G6OsYmJjyISlaN4+B056Ctjd6eXq+jqUkT4xNwh4kx22CL94uIzCM3Pk5rS8uc443G4CsUmCgUlJzVumSyfDX87CCZoJ9I+PaMqRKzpcuP5zENZs7JD4C6QB1j18Y8iGpplJzVkI0tLWT8fobzedr9M//Tnrp1i65t2zyKTFZNIEAooAWNKyW0KQQfFpdPTB8ArLWYa4b2h9vv8mwRWcu2b9rEh1eu8PCssvCDY2M0t7bS5NdPspqWTBL60gCU2/b46Y0+72lItSS4LojJGibHJ6lrmHmSY+zyGJu2VMesGSg5qykNdXU8sW8ff/jOO/xwezvBujqstXyUy/FqXR2xcNjrEGUFxFOJ4pUQpIJZQtpttmLWd61nXcM6rvdfp21vG6bOYCctwx8Os6FpA+u2rfM6RBGpUEf27uWPL15k/dgYWxuLU/BXJyb4k1u3ePqRR+Y94y+1Ix7JaGXLCmpoamD3/t30v9VPe6Sdusbib+DRC6M0XGlg+2eqZ8uBkrMa85kHH6RgLb988iSbJicZtZZCays/9OijbJpnOYVUsVJlp9TBPJTOxIYCXdprtoKMz/Dkdz3JOy++w4UXL0ALcBO6tnbxyHfqx5WI3NkDHR184amn+P2336Z5eBg/cK2+nk8/9hiHtm71OjyRqtf9ZDc+4+PUsVMUWgrYMUuoOcSj3/0ogWD1VENVclZjfMbwuT17+NSuXVy8cYNGv58t82xAlirlOMRDqeLV7ly5spPWqa+exuZGHv/Ox7l14xa5GzkCrQGaWpu8DktEqsDBzZvZ//nPc+HGDQrWsrW1VfvM1oBysS7avA6lpvnqfHQ/1c2Dhx4kez2Lv8FfXO5YZb+BlZzVqIDfz66ODq/DkOVUquyUCfohEAACRDojSsw80tTapKRMRBatzudje7v2p64VU4lZaNMOrWxZJfWBejq2VO9vYCVnIpUsmSyuUweSz2ZLlZ16lJCJiIhUiza1t5GFU3ImUqmmKju1TS2DCGojsYiIiEgNU3ImUkHiqQROZw6A1JfGVdlJREREZA1RcibiNccBIB5KFTcMt7URChQL4msZhIiISJVyHJzuHFA9lQLFe0rORDwUTR+FSPFLOxXMEtqgDcMiIiJVb6rdTSdEOiNeRyNVRMmZyGoqzZIBxS/tDXlCG4pNo0OElJiJiIjUgHJipiJeskhKzkRWS+ksmlvMxcj41aNMRESkVoVCnRrjZdGUnImspKmZMtct9yiLhHvKd+tLW0RERESmKDkTWSmlptH4/dADmYB6lImIiIjInSk5E1lOyWT56lSPsqmNwErK7i0/nmfgvQHOnjzLxPgEm7ZsYndkN20b2+79ZBEREVkR1lounrrIR+99RHYkS1uojd0P72bTrk0YY7wOr6YoORNZJvFUguSzWQgGyZCDgHqULcbkxCR9/7OPK/YKzfua8Qf8DF4cZPBPB/nUt3+K9dvXex2iiIjImtT/Sj8fnP2Apj1N1HfXM5wZ5pVjr3Dw2kH2HN7jdXg1RcmZyP0ozZTFI5lyjzLNlC3NhQ8vcCV/hdDhUPksXNsDbdxqvcVbf/cWn/vi53R2TkREKl48lSB1YJyQ14Esk+y1LCc/PEnoSAhfvQ8Af5OfxnWNvH/sfbr2ddHU2uRxlLVDyZnIEkXTR0n9o3xxTxkQCqlH2f0YOD1AYGdgTgIW2BBg+INhsleztG5o9Sg6ERGRe4unEuWTtbXym+DS2UvQSTkxm1LXWIfdaBk6N8T2g9s9iq72KDkTWSjHIR5KFa925ko9yrpq5svXa5OTk5i6uTNjxhhMnaEwWfAgKhERkYVzOnOENtXWydrJyUmou8OddWh8XmZKzkQWwnEIHz5GJuiHQAAIEOmMaOniMtqyfQvvXXiPwPrAjOMT2QnqJ+oJrg96FJmIiMgCOA5EAvd8WLXZsG0D9qTF7rEzVrfYgoUr0PFEh4fR1Z5VTc4ujF4mnkrMOBaPqGCCVKhkkngkU7x6OFvuUaaEbGXsOLCDM/1nGPlohNZwK6bOMD48TvadLI8+9ih1/judthMREZGV0rG1gy2hLVx8+yJt+9uoa6wjfyvPjfdvsGPLDm05WGarmpxdDkKiZ9qBkRFIJZSgSeUp9Si7PVMW1EzZCmtsaeTpv/c0J14+wcW/vQh+aPI38fhjj9O1v8vr8ERERNYkYwyHv3CYk986yZljZyj4C9QV6tjfvZ89j+9Rsa5ltqrJ2aaWTTzX81z5djKVJOEbxEkfJZq+v2ngeCoEsdh9Rihr2fRZ3cSXipt5VQp/dTW3N/PEdz3BRG6C/ESeQEsA49OXvoiIVDjHIdrdRyqYJ1QzdRpv8zf4eejIQ3Q/2c34rXEamhqoq9eKlpXg6Z6zWCRGMpUkFcyQCt/fayUODJBJJpWgyZJMr65UpMTMS/WBeuoD9V6HISIicm9TiVknRMJHanqVTV19HU31Kpu/kjwvCLJc1WwSfQlCXxogkj5636/l9PdANHr/QUlFi6cSOJ05gGI/khqrriQiIiKrJBAgFArVdGImq8Pz5Gy59Pb0kkwlcTfd3+tkMmnCwWO4zrSDStRqg+OUr0a7+0gdLJbCBwixfCcKRERERESWomaSM1ieH9eO65By+wj39BHOBSCXw0m6Wi5Z7UpLDorFPSiuCVePMhERERGpIDWVnC2HqeloN+PiUppJe3YQV/vZqs+0mbKpHmWhUHGTbiSkyosiIiKyDFwXtyfndRRSI5SczWP2j/ap/WwqOFJFSqXw8Rf/iWcC6lEmIiIiyyyZJPSlgWKFZ63GkWXg8zqAatDb0wttbYS+NDCnibZUEMcpXqb1KAt39xDu7lFiJiIiIsvLcYg+k1brHVlWmjlboN6eXhJ9pXLrapxdceKpBMnDWfD7yUSBQFBflCIiIrKyAgFCgdrraybeUXK2CFMVIRMMKEGrBMkkAPFIptyjLNIZAeYuTRURERERqXRKzhZpqnF2wjeIkz6K0/m81yGtSfFUguSzWQgGyZAjFFKPMhERERGpbvfcc2aM+S1jzGVjzIlpx+LGmI+NManS5btWNszKEovECG3oIrUhTzR9dEZVQFkhpb1kJJNE00dJHMxCVxfhcIRIuEeJmYisSRqjRURqy0JmzpLALwO/Pev4L1prf27ZI6oSsUis2BONY0S7+4r5mZpVr4xSj7LUkXy5+mIopB5lIiJojBYRqSn3TM6stS8ZY8IrH0r1mdrXlPL3EQ4ew3VQgrbcHKfcoywSPqK9ZCIi02iMFvHI1InjYJ4QKggiC5fou3vl9/sppf+Txph3SksqOu7jdapaNBwlEu4plm0/fKxcpEKWQTI5LTFTKXwRkUXQGC2yUqYSsw15Qhu0kkcWLtGXgJGRuz5mqQVBfhX494At/e/PA/90vgcaY54DngNo39y+xLerbOUZNLeP8LODuGpWvTzCYWDQ6yhERKrNgsbo6ePzjvbaHJ9FVkwgQGhDSImZzMtxHVLp1JzjoaEs7gs7MJy743OXNHNmrb1krZ201haA3wCeuMtjv2qtPWytPdzc3ryUt6sK0XCU3ujzZDYECT87qBm05RCN4h4/QiQNqdPHcFzH64hERCreQsfo6ePzxubaHZ9FRFaT4zqk3D5CQ1l6+5hxcV/ouucEzpJmzowxW6y1F0s3fwA4cbfHryWRzghubm6mLEsUjeI4FJcPcAw34+oslYjIXWiMFhFZXVMJGQD5PJEhP07/kbm1KCL3fq17JmfGmN8DosAGY8wg8DNA1BgTobhkwgV+YuHh174MOeKRDHGvA6kVMxK0QRzX0f4zERE0Rot4wnVxe3JeRyEem1rR5WZcMkODxYTsG53FO8PhJRcJXEi1xi/Oc/g3l/Rua0A0HCWVTpE4MAKpBPFIr9ch1YZSghbu6fM6EhGRiqExWmSVJZOEnx0kE/DTq5U8a1YylSSTSRMiALkckYwfp/N5iN3/ay+1IIjcRW9PL8lUkgQDStCWWz6Pm3G9jmLJrLVcv3CdywOXAdi0YxMdWzswxngcmYiIyNp2KZvlrcsXyebH2d22noc2bqS+ru72AxyH6DNpMhuCRDojWsmzxkzNlKXSKRgZoff9NuKpUhuFZSwEqORshcQisWKC5hvESR8tZtNyf6JRYqkUiYYBEn0JenuqK+mdzE/y5jff5OOhjzGdxWSs/2/62bZhG499x2P46u6ns4WIiIgs1V+ePc3vXzyB2WKoazH82eVTdLlt/NQnnqKjqan4oGiUaCqFG8ziZvrI+POk0qmq+z0ii5foS0A2SyjvJ5TPEzvdVpx8WcAessVScraCphK0FINE00dx+nvUpPo+xSO9kEqQOJglmUpWVXGQs6mznL95no5PdWB8xeTMPmg5/8Z51qXW8eBjD3ocoYiIyNrz0bVr/N7lE2x7vJX6+tJM2Tb4+PwIyf4U/+rQU+XHxiO9xB2neD1U3MaS6EsQ6Zz5K10zatVpvsrgUzNlma/vKLV5YkV/zys5W2EzErTuPhwHJWj3KZ6J4GRTuBu8jmThrLWcfu80rY+1lhMzAOMzBPcFOZ06reRMRETEAy9ePEfD9rrbiVnJ1q5W3jl/iaHRUTZMbzdR+h0XJ0o8mST6TBo3M3NPfMLtIxLuUZJWRcol8GfVeolkwPnGjlXrYazkbBXEIrHif3COKUFbLrkcmVzG6ygWzBYsuVyOUDA057761nqGR4exBTsjcRMREZGVlx7L0tJSP+e4MYa6ZsPI2NjM5Gy6WAynNJNW5rqEnx0k5faV98mHQ2ElahXIcZ3yf6NMJk0om8c9fmTuA2PRVYtJydkqmfo/ZMrfRzh4DNdBCdpSRaM4yeIXX8I5Sm+08vfzGZ+hpbWFsetjNHY0zrhv/Po4re2tSsxEREQ8sKs5xJmR67S3B2Ycn5wsULgJ66f2nN3JPL/nXMch2t0HmeKJ5FTwWPGhStAqRjKVLJbAzwaLB3LM35tslSk5W0XlBM3tI3z4GG7SXbUp0poTiq9DOAAAHhhJREFUi+Emk4T/Yboq9p4ZY9j3yD6Ov3Oc+sfr8dUXi38UJgrc/OAmjx963OMIRURE1qbPbAvzzXfPMLp+gubm4gyatZbzZ0Z4un0H7YHAPV5hHqUWQOWb3X2kODaj4nSl/3aZbmqGqZpink8ylSxfL/cm65+2X7ACJk6UnK2yaDhKNBwl4Rwl/OwgbjKpBG2pwmHCuQyu13Es0PaHtnNz5CYnXzoJpf1yZshw4MABuvZ3eRuciIjIGrW1tZX//YEn+PU33+BKx01oAHsNDjV28qWHHl76C0/7oe8QJZo+CkNpAFLrxqvi5DJM24uVzZPIVMeKpfkkU0kylweIXGsAIJoOFgvNdXoc2CzGWrtqb7Z131b73K8/t2rvV+kSfYnb1V+UoC1eaclAqpOq2nR7a+QWVwevYoxhXdc6mlrvsVxCKprjOsVKThVo+N8Ov2GtPex1HFL5Dm/dao8/p/FZ1rZbExP0Dw2Ry+fpamujq61txfqQ/v/t3W1sXNd95/HvIYfi04gaihJFPdG0ZNmyYytj+SGyrCjjuMmmybbpInG6CVJk2kXdAl0g3L6pdt+UKVBAWHR3580iRbLNTnYRtwmUJummhZEg9oWtpHQs2zeyJTMSbY0pShpJjDSihuKIHPLsi5mh+CiNyCHvvcPfBzBI3qHJ/8EV58xvzlOPmyDx0DC0tExd89OW/DP6tlyOaBqcFzsKB3CHQ7CY0USvlc4m88H5w+ZrX1uwf9bImYe693WT6E0Q+dIA3Tqs+u4VpwzE6MXNB2cud2NLI9se0khZNZj+bmK8P+x1OXN8zesCREQCpLGujkc3b16R31U6GqiktCW/HwLa1FqsoRCxdAPQQE8mCvEYqWSSnmhwNmSbyR/B7E4UzjxWCmiySKWANm0udxCmCEjwJN3k3B1Cc7lbOzv5YJ76bF/7keKZiIhfTQ8KPckkkS8NzPuaMNoRXfY3n2f0cdlsYS1Wx6G5U/7icXqWtRJROJPgmxHQBgMzh1uCY+67iCXFdxN9GMxERCRA4nEy84xKJe/L4uYKZ6gtV0ArLbPpPlmaYhkOxAhTtVI4k+pQCmjRYB1OLf5V9ruIIiIilTDPqFSP49D1+FHcVO/UGrBKjKTNGKHz0VosUTiTapPL3fl7RO5g9ruIPe42bdojIiIrLxYj5UBPpBDMnI4cbnZpyzjm9nHamM5PFM6kesRixFyXxPoB3yyqleBI9CZuhfuxsZnvIkYX/v9ERESWVSxGD7FbX6YP4zJIwjkMQCTScdugVtq8qiSSzZM6Mi2QqY/zFYUzqSql3Y8SDw1r7ZnckZNyAKZ2XEwdmbaLpt5FFBERH3I6DkEyCUBPNEPioQGSbpKuSNe8368+LlgUzjyWcA4TyeYL06b0zkVF9ES7cdKHSbV7XYn41dQW+BQ295jqtNRZiYhIEBT7qx4AN0FybJBUZp4t7nM5ohlw+g5APLZy9cmiKZx5KNGb0IvCZTRn23NZstJI0+347ay5+Wq+9S5ipHglor9BEREJpJ5oNz3FkbR5dXVpV+EAUTjzWLw/rBeFy6B0in3COUx37JDX5VSFRG8Cslki+YWfNjKhPG7a9c16v9IW+LNr1ruIIiJSVfRasmr4JpxZa7ly7gpnf32Wmzdv0rapje27t1PfXO91aRJE8TipZLIQ0LQ5yKI5KYdUJlUYhRweJvNCZ+EduIWkUlOhOBLpoCvSteIjaaWagamzyZy+fXO/Ue8iipTtYjbLz8+f5ezoNbY2ruXpzZ1sXrvW67JERKqOL8KZtZaTR09y6swpQp0hatfXcuH8BU6/fZqnP/M0LRtb7vxDRGaLx4m7CRI69+yuTa3JykFXrgFyOZwXy9tqN5VMEvtUGjIZUg29JFK9RLv2rUhIK43uRbNhyOWIpYsHaepsMpFFezN9gf/53i+xWyG8YQ3Hr1/kX97u50+6HuMjW7bd+QeIiEjZfBHOfnP2N5xKnWLd/nXUhGoAaOpoYuTcCG+8/Aax52IYYzyuUgIrl8NJOb5bC+U3M0acMulb6yFLI2XlTgGMx3EcZ+rL0uGZpZ9d6dG0pJsEmBrd6z7ZQk+muLuORsdElmRkbIy/7T9G294mmprqAFi/vpHR9nH+1xtv8mDbRlrqNcNFRKRSfBHOBn49QKgzNBXMSpq2NHHtvWtcH7qu0TNZlJ5MFCfdi5s/Cvhvswq/cFIObv/RwogTQG6Ja7KmhaKUA7HdvVDcRcoNL+3wzOlmj5SVO7onIuV559IlbrZNTAWzksbGOsY2Zjl+Mc2Bzns8qk5EpPr4IpzlcjlCkbmlGGMwDYbxm+MeVCVVIRbDcQrhwKVyoaAalEacYPrarGkjTpWaChiL4RCD4mha4V4Mzvj95d6Tio3uiUhZbuTHMQ123sdqGgzZcfXPIiKV5Itw1tbexuXLl2nc2Djj+uT4JAxDuDXsUWVSFYrhIJY+PBUKVntAS7pJMpcGiF5ZU7wSKhxquZxrs4qjaQ4xetwETkcaAHf9WFn3ZGodXDZPV7bw1OW8qGMoRJbTlvBaGDBw79zH7FXYtkWzWkREKskX4azzwU76v9/P6IZRGtsLAW0yP8m149fYed9O7dgoFeF0HKLHTZDYtzrPP0v0Jm59UVqbFfVmF8vpv7fHTZBgYGZ988lmi6N7B25Nm4wvW4kiAuxqa2PH+6188EGGbZ0tGGOw1nJ+8Drbb67jwQ3acUlEpJJ8Ec4aWxrZ/6n9vPHSG2ROZzD1Bq7Bjh07+NCBD3ldnkjgJXoTU4GswLtgNltPtBvcOwQzAMKFjT60yYfIiqkxhq/u+Qh/9+5bvH3hIjVhw2TW8mDdRv54z15qa2ru/ENERKRsvghnAK2bW3n2i89y7dI1xm+Os7ZtLQ3hBq/LWjalF8ugKSErbnh4VZx9NnukLPOCfzfL8EtQFJG51jU08OePPkU6m+XK6Cit9zTQEQ5rF2URkWXgm3AGYGoMkY6I12Usu+mjGHpRurJKozSJh6o3oE1fmxXvL6zX7HH9G8xEJBg6wmE6wloDLiKynHwVzlYTBTPvTAW0fV5XUjlOysFNu4UvcrnCLobHpq3NinpWmoiIiIiUSeFMJOCmzigbCuG8WNxusatLa7NEREREAkbhTFalHjdC8r5BEs5humOHvC7nriXdJJlMYSt68vniLob7dM6XiIiISIApnMnqFI+TSibp+vxgYNaeOSkHgFQmRebSQGFqrFtcoxmPL+8ZZSIiIiKy7BTOZPWKx4m7CRIBOKYn0ZsorCWjAbJZuvuLaxa1lkxERESkaiicieRyOCmHWFfM60pmKI2UuWn31u6ebgSIaOdFERERkSqkcCarWk8mipPuxc0fBfBNQEs4h4nkCp9H8nlSR4pb4WukTERERKRqKZytMCflQDYL6KwYX4jFcByI7e7F5SipTIp4NL7iZTgph1QmBUAmky5shX9kW2HXRdBGHyIiIiKrgMLZCkq6STJDg0SHQvRkNATiG7EYDjFi6cO4DJJ0kysa0KYOjc5BV64BcuD0HVAgExEREVll7hjOjDHfAv4tcMla+3Dx2nrgu0AXkAK+YK29unxlVodMLnNry/MqPYPqYjbLv545w+ClS9SvWcOeri4e27qVUE2N16XdkdNxiB43QSKaXvY1aEk3OfV5KbA7fdNOxa7Sfx8iUlnqo6VcI2Nj9A4McGpwEGst92/bxkc6O1lbX+91aSIyTTmvmJPAp2ZdOwT8zFq7C/hZ8WspQyzdULUvvN+/epVvv/wyG95/ny9MTPBsNsvpY8f4zuuvk5+c9Lq8svS4kcKOiMuoNILalcrQlZoV2Ev/iYiUJ4n6aLmD4Zs3+earrzL6zjv87s2b/N7YGOMnTvDNV14hk8t5XZ6ITHPHkTNr7SvGmK5Zlz8LxIqffxtwgL+oYF0SMNZafvzmm3wuFGJnYyMA7cC9DQ383/Pn+VU6zWNbtnhbZLmyWdy0W9GRsxkjZaUzykpTW2MxnVEmIouiPlrK4fT388j16zwbiUxd21JfT9O1a/zs17/mcx/+sIfVich0i51rtslaewGg+LG9ciVJEF3IZqnNZtnRMHPUyRjDk/X1nBgY8KiyuxSPE+8Pw/Bw4WyxCkj0JshcGqCrL01XX7oQzKLdGiUTkeWiPlqmWGt558wZPrJ27ZzHngyHefeDD5i01oPKRGQ+y74hiDHmeeB5gHWb1i33rxOPjE1M0EghjM3WVFPD2Pj4yhe1SD3RbnqSSSLxxa89mxHsSmeURbsLX2uUTER8YHr/3LlO/XM1G8/naZxn7Xd9TQ12cpJJa6mZp/8WkZW32JGzi8aYzQDFj5cW+kZr7TestY9bax9vWte0yF8nfrc5HOZSbS3X8/k5j72by3FPR/ASyWLXniWcw0SGsnT3QncvZF7ovBXMRESWX1l99PT+eWOT+udqZYyhc9Mm+m7cmPPY6dFROjZsCMSmXSKrxWJHzv4J+ApwuPjxRxWrSAKpPhTiyd27+d7x43yupYVIKMSktRwfGeF4fT3Pb9/udYl3L5vFTfUCtz+c2kk5uGm38EUud+uMsni8cE2nJojIylIfLTMcfOABfvjKK4RzOTrr6zHGcDaX45/HxvjM7t1elyci05Szlf7fU1hYvMEYMwj8JYUn/O8ZY/4DMAA8t5xFSjA8s3MntTU1fKOvj7UjI4wAkQ0b+IMPf5h1Dcu7A2LFxeOkkkm6Pj9424A2dUZZNl9Yq0ZDYaMPnVEmIitAfbSUY+f69Xx6/35+9PbbMDyMASaamvjkU0/xwIYNXpcnItOUs1vjFxd46NkK11LVEs5hItk8Pe62qh1JMcbwsR072H/PPQzduEFDKERrcefGQCoFtHhmxuWkmySTSRe+yOeLW+Ef0OYeIrLi1EdLuR5sb2f3xz/O5Rs3sNaysblZ68xEfGjZNwSRwuYQc6a6VbG62lo2z7MrVDUonVEWHQrhvFhcR9fVpWAmIiK+Z4yhvbnZ6zJE5DYUzlZCLrdqglnV6eqiK92Lmz9amN44NnZr58W418WJiIiISDVROFtmTsrxugRZilgMxwFSqVvXFLJFREREZBkonC2j0mYR0TSFqW8STJqyKCIiIiIrQOFsmcxYm9S3Ty/wRURERETktnTq4DJQMBMRERERkbulkbMKmxHMOg5Bh9cViYiIiIhIEGjkrIKSbpLMpQG63wkXgpmIiIiIiEiZFM4qxEk5ZDLpW9usi4iIiIiI3AWFswqK0ECPG/G6DBERERERCSCtOasQN+0SyWYBhTPxp3Q2yyunT9N/7hzGGHZv387Hdu1ifWOj16UFyqS1vD44yOunT3M1m6U1HObJ++/nia1bMcZ4XZ6IiATMyNgYr7z/Pm+fOcPY+Didmzbx0fvv597WVq9LC5zTv/kNR0+d4tzlyzSsWcOeHTs4uGMHDaHgRJ7gVOpjid4EkaEsqSPbdECx+NK54WG+8+qrfGxigt9pbmYCePPMGb514QJ/dPCgAtpd+H8nTnDl1Ck+29TE5rVrOX/zJj/95S+5+MAD/M6HPuR1eSIiEiCj4+N86xe/YOfVq/xxOExzfT19ly/z/XSazzz1FA+2t3tdYmC4Fy7w8muv8W/q6rgvHOb6xARHT5wgefEif7R/P2tqa70usSya1rgETsoh4RxWMBPf+9m77/KJyUk+0tJCY20t4dpaDq5bxxO5HK/093tdXmBcuH6d9/r7+XIkwvaGBkLG0NnQwB9EIpw6fZpLIyNelygiIgHy+uAg265e5dOtrbTW1bGmpoY94TDP1dfzk+PHmbTW6xIDIT85yU9/9Su+1NTEQ83NrKmpoa2ujt+NRFg3NIR74YLXJZZN4WyRnJSDm+olks0rmImvjU1McDad5pFweM5je5ub6Tt71oOqgunXQ0PssZa6mplPnWtqanjYWvouX/aoMhERCaK+s2fZO8/slc76empGRvSmX5kGh4dpvXmTTWvWzLhujGFvfT19g4MeVXb3NK1xEWYEs2MHIB7zuiSRBVlrwdp534mpNYbJyckVrymoJq2lZoF1ZbXFx0VERMplraV2nn7FGEOttepXyjS5wOscCN5rHY2c3SUn5eD2HyWaphDMYjGvSxK5rfpQiI6NG+m7cWPOY8dHRti1dasHVQXTfW1tnAAmZnWWE9Zywhh2tbV5U5iIiATSfVu3cnx0dM719NgYo42NbGpu9qCq4NnW0sLlUIir4+NzHjueywXqtY7C2V1IuslCMBsK4fTtUzCTwHjmoYf453yekyMjTFpL3lrevH6dV2pr+eiuXV6XFxjbW1rYsH07RzKZqQ7gyvg437t6lY577mFrS4vHFYqISJA8uX07fc3NvHrtGjcnJ7HWcmZ0lO+OjPDMI49QW6OX6uVYU1vLwUce4YXr1zmby2GtZXRigpeuXePsunXs3bLF6xLLpmmNZUq6STJDg4Vg1nEIOryuSKR8O1pb+dzBg7x88iQ/uHwZjKFz82a+/OCDdMyzFk3mZ4zhuUcfxWlp4Zv9/UyOjFBTX89je/bwsR07vC5PREQCZm19PX944AA/7evjb86epWZyknWRCM/u3cvDmzZ5XV6gPNXZSUMoxA/7+sheu4atrWX3vffyhw88QGNdndfllU3h7C5Es2GcvqiCmQTSjtZWdjz9NDfzeYwxgdlS1m9CNTX81q5dfPy++8jl8zSEQguuQxMREbmT1sZGvvDoo4zv2UN+cpKGUEjnZi7So1u2EN28mVw+T11tLaEAjjwqnN2NXM7rCkSWrD5ABzH6WY0xNAXonTgREfG3utpa6vTG6ZIZYwI1UjabXqWVIdGbgOFhYukWrTMTEREREZFlEbyxvhVWCmbdJ1voiXZ7XY6IiIiIiFQpjZzdRimYZV7o1CHTIiIiIiKyrBTOFpBwDhcOmT6iYCYiIiIiIstP4WwWJ+XgpnqLwWybgpmIiIiIiKwIhbNZUplUIZgdOwDxmNfliIiIiIjIKqENQebRlQ1pV0YREREREVlRCmfTOCmHzNCg12WIiIiIiMgqpGmNRUk3SWZokOhQCKdvH3R4XZGIiIiIiKwmCmfMCmYdhxTMRERERERkxa36aY1JN0nm0gDd74QLwUxERERERMQDq3rkrHTIdPfJFnqi3V6XIyIiIiIiq9iqHTlTMBMRERERET9ZlSNnpWCWeaFTh0yLiIiIiIgvrLqRs0RvgshQVsFMRERERER8ZdWFM4B4f1jBTEREREREfGVVhTMn5UAu53UZIiIiIiIic6yacOakHNxUL5Fsnp5M1OtyREREREREZlgVG4I4KQe3/2jhkOm+AxCLeV2SL1hr6b9yhZPnz5OfmKCrvZ09mzZRV1vrdWkiIiKr2sVslrfOneP6jRu0RyLs3bKFtfX1XpclIstsSeHMGJMCrgMTQN5a+3gliqqkpJskMzRYDGb7FMyKJq3liOvymw8+YG9NDfU1NZxMpfhFJEL8qafUAYiIBFwQ+miZ378ODPDzt97icWBrbS0DAwN8/d13+cL+/XRFIl6XJyLLqBIjZ89Ya4cq8HMqTsFsYW+cO8eNM2d4vrWVWmMAiAIvXbvGv5w4we/v3ettgSIiUgm+7aNlfpdGRvj5W2/xfHMzLaHCy7RHgAdHRzny2mv8p098gtqaVbMqRWTVqdq/7lIw634njNNxSMFslrfef5+DjY1Twazk6bVrOTM4yI3xcY8qExERWb3cc+fYC1PBrGRHYyNto6OcvnLFm8JEZEUsNZxZ4CfGmDeMMc/P9w3GmOeNMceMMcduXLuxxF9XHiflkMmk6X4nTE+0e0V+Z9BkR0dpq6ubc72+poYmaxlVOBMRCbrb9tHT++fLN1amf5Y7y46O0rbA2u82YGRsbGULEpEVtdRw9rS1di/w28CfGWMOzv4Ga+03rLWPW2sfb1rXtMRfV74IDfS4mpe9kPbWVlLzHCuQyecZDYVo0ZozEZGgu20fPb1/3ti0cv2z3F77unWk8vk51621fAC0NzevfFEismKWFM6steeLHy8BPwCerERRS+WmXchmvS7D15667z5eGh/nyrQRsrHJSX58/TqP3X+/dmwUEQk4v/bRcnuPbtnCqfp6Tk0bzbTW4gwP09jezraWFg+rE5HltugNQYwxzUCNtfZ68fNPAn9VscoWKdGbgOFhUkc6IR73uhzf2rl+PR994gm+4bp0joxQD7xnDLvvv59ndu70ujwREVkCv/bRcmfNa9bw7/fv58jrr/NqJkObMQxYS3N7O7//2GOYWWvFRaS6LGW3xk3AD4pPEiHgBWvtixWpapESzmEi2byCWZke37aNRzo6OH3lCvnJSZ6NRIg0NHhdloiILJ3v+mgp3/Z16/jqs8/y3pUrZMfGeKK5mS1r1yqYiawCiw5n1tr3gQ9XsJZFc1IObqq3GMy2KZjdhfpQiIfb270uQ0REKshPfbQsTo0x7Gpr87oMEVlhgd9KX8FMRERERESqQSUOofbMjGB27ADEY16XJCIiIiIisiiBHTlzUg5u/1GiaQrBTIdMi4iIiIhIgAUynJVGzKJDIZy+fQpmIiIiIiISeMZau3K/zJjLwAdlfOsGYGiZy1lJ1dSeamoLqD1+Vk1tAW/ac4+1duMK/04JoLvon6G6/jarqS1QXe2pprZAdbWnmtoCPuufVzSclcsYc8xa+7jXdVRKNbWnmtoCao+fVVNboPraI6tXNf1brqa2QHW1p5raAtXVnmpqC/ivPYGc1igiIiIiIlJtFM5ERERERER8wK/h7BteF1Bh1dSeamoLqD1+Vk1tgeprj6xe1fRvuZraAtXVnmpqC1RXe6qpLeCz9vhyzZmIiIiIiMhq49eRMxERERERkVXFd+HMGJMyxrxtjHGNMce8ruduGWO+ZYy5ZIx5Z9q19caYnxpjThc/tnpZY7kWaEuPMeZc8f64xphPe1ljuYwx240xLxtj3jXGnDDGfLV4Paj3ZqH2BPX+NBhjfmmM+VWxPV8rXr/XGPNa8f581xizxuta7+Q2bUkaY85MuzdRr2sVuRvqn/2jmvpnqK4+Wv2zfwWlf/bdtEZjTAp43FobyPMTjDEHgSzwf6y1Dxev/VfgirX2sDHmENBqrf0LL+ssxwJt6QGy1tq/8bK2u2WM2Qxstta+aYxZC7wB/B4QJ5j3ZqH2fIFg3h8DNFtrs8aYOuAo8FXgz4F/tNb+gzHmb4FfWWu/7mWtd3Kbtvwp8GNr7RFPCxRZJPXP/lFN/TNUVx+t/tm/gtI/+27kLOista8AV2Zd/izw7eLn36bwR+p7C7QlkKy1F6y1bxY/vw68C2wluPdmofYEki3IFr+sK/5ngY8DpSfLQNyf27RFRDyk/tm/qqmPVv/sX0Hpn/0YzizwE2PMG8aY570upkI2WWsvQOGPFmj3uJ6l+o/GmOPFaRW+n2IwmzGmC3gUeI0quDez2gMBvT/GmFpjjAtcAn4KvAdkrLX54rcMEpAObnZbrLWle/PXxXvzP4wx9R6WKLIY6p/9L5DP/9NVUx+t/tl/gtA/+zGcPW2t3Qv8NvBnxaF78Y+vAzuBKHAB+G/elnN3jDFh4PtAt7V22Ot6lmqe9gT2/lhrJ6y1UWAb8CTw4HzftrJVLc7sthhjHgb+M7AbeAJYD/h6ao7IPNQ/+1tgn/9LqqmPVv/sT0Hon30Xzqy154sfLwE/oPCPIOguFucgl+YiX/K4nkWz1l4s/sOeBL5JgO5PcX7x94HvWGv/sXg5sPdmvvYE+f6UWGszgAPsAyLGmFDxoW3Aea/qWoxpbflUcaqLtdbeBP43Abw3srqpf/a3oD//V1Mfrf7Z//zcP/sqnBljmouLJzHGNAOfBN65/f8VCP8EfKX4+VeAH3lYy5KUniSL/h0BuT/FRaB/B7xrrf3v0x4K5L1ZqD0Bvj8bjTGR4ueNwG9RmKf/MvD54rcF4v4s0Ja+aS8wDIW5+YG4NyKg/jkIgvr8D9XVR6t/9q+g9M++2q3RGLODwrtxACHgBWvtX3tY0l0zxvw9EAM2ABeBvwR+CHwP6AQGgOestb5fyLtAW2IUhuQtkAL+pDQf3M+MMQeAV4G3gcni5f9CYR54EO/NQu35IsG8P3soLCiupfCm0festX9VfE74BwrTDN4Cvlx8Z8u3btOWl4CNgAFc4E+nLUwW8TX1z/5STf0zVFcfrf7Zv4LSP/sqnImIiIiIiKxWvprWKCIiIiIislopnImIiIiIiPiAwpmIiIiIiIgPKJyJiIiIiIj4gMKZiIiIiIiIDyiciYiIiIiI+IDCmYiIiIiIiA8onImIiIiIiPjA/weyvjipu20WGQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1080x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "titles = ('K Neighbors with k=1', 'K Neighbors with k=2')\n",
    "\n",
    "fig = plt.figure(figsize=(15, 5))\n",
    "plt.subplots_adjust(wspace=0.4, hspace=0.4)\n",
    "\n",
    "X0, X1 = X_train[:, 0], X_train[:, 1]\n",
    "\n",
    "x_min, x_max = X0.min() - 1, X0.max() + 1\n",
    "y_min, y_max = X1.min() - 1, X1.max() + 1\n",
    "xx, yy = np.meshgrid(np.arange(x_min, x_max, 0.2),\n",
    "                     np.arange(y_min, y_max, 0.2))\n",
    "\n",
    "for clf, title, ax in zip(models, titles, fig.subplots(1, 2).flatten()):\n",
    "    Z = clf.predict(np.c_[xx.ravel(), yy.ravel()])\n",
    "    Z = Z.reshape(xx.shape)\n",
    "    colors = ('red', 'green', 'lightgreen', 'gray', 'cyan')\n",
    "    cmap = ListedColormap(colors[:len(np.unique(Z))])\n",
    "    ax.contourf(xx, yy, Z, cmap=cmap, alpha=0.5)\n",
    "    ax.scatter(X0, X1, c=y_train, s=50, edgecolors='k', cmap=cmap, alpha=0.5)\n",
    "    ax.set_title(title)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 习题3.2\n",
    "&emsp;&emsp;利用例题3.2构造的$kd$树求点$x=(3,4.5)^T$的最近邻点。\n",
    "\n",
    "**解答：**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x点的最近邻点是(2, 3)\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "from sklearn.neighbors import KDTree\n",
    "\n",
    "train_data = np.array([(2, 3), (5, 4), (9, 6), (4, 7), (8, 1), (7, 2)])\n",
    "tree = KDTree(train_data, leaf_size=2)\n",
    "dist, ind = tree.query(np.array([(3, 4.5)]), k=1)\n",
    "x1 = train_data[ind[0]][0][0]\n",
    "x2 = train_data[ind[0]][0][1]\n",
    "\n",
    "print(\"x点的最近邻点是({0}, {1})\".format(x1, x2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 习题3.3\n",
    "&emsp;&emsp;参照算法3.3，写出输出为$x$的$k$近邻的算法。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**解答：**  \n",
    "**算法：用kd树的$k$近邻搜索**  \n",
    "输入：已构造的kd树；目标点$x$；    \n",
    "输出：$x$的最近邻    \n",
    "1. 在$kd$树中找出包含目标点$x$的叶结点：从根结点出发，递归地向下访问树。若目标点$x$当前维的坐标小于切分点的坐标，则移动到左子结点，否则移动到右子结点，直到子结点为叶结点为止；  \n",
    "2. 如果“当前$k$近邻点集”元素数量小于$k$或者叶节点距离小于“当前$k$近邻点集”中最远点距离，那么将叶节点插入“当前k近邻点集”；  \n",
    "3. 递归地向上回退，在每个结点进行以下操作：  \n",
    "(a)如果“当前$k$近邻点集”元素数量小于$k$或者当前节点距离小于“当前$k$近邻点集”中最远点距离，那么将该节点插入“当前$k$近邻点集”。  \n",
    "(b)检查另一子结点对应的区域是否与以目标点为球心、以目标点与于“当前$k$近邻点集”中最远点间的距离为半径的超球体相交。如果相交，可能在另一个子结点对应的区域内存在距目标点更近的点，移动到另一个子结点，接着，递归地进行最近邻搜索；如果不相交，向上回退；\n",
    "4. 当回退到根结点时，搜索结束，最后的“当前$k$近邻点集”即为$x$的最近邻点。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 构建kd树，搜索待预测点所属区域\n",
    "from collections import namedtuple\n",
    "import numpy as np\n",
    "\n",
    "\n",
    "# 建立节点类\n",
    "class Node(namedtuple(\"Node\", \"location left_child right_child\")):\n",
    "    def __repr__(self):\n",
    "        return str(tuple(self))\n",
    "\n",
    "\n",
    "# kd tree类\n",
    "class KdTree():\n",
    "    def __init__(self, k=1):\n",
    "        self.k = k\n",
    "        self.kdtree = None\n",
    "\n",
    "    # 构建kd tree\n",
    "    def _fit(self, X, depth=0):\n",
    "        try:\n",
    "            k = self.k\n",
    "        except IndexError as e:\n",
    "            return None\n",
    "        # 这里可以展开，通过方差选择axis\n",
    "        axis = depth % k\n",
    "        X = X[X[:, axis].argsort()]\n",
    "        median = X.shape[0] // 2\n",
    "        try:\n",
    "            X[median]\n",
    "        except IndexError:\n",
    "            return None\n",
    "        return Node(location=X[median],\n",
    "                    left_child=self._fit(X[:median], depth + 1),\n",
    "                    right_child=self._fit(X[median + 1:], depth + 1))\n",
    "\n",
    "    def _search(self, point, tree=None, depth=0, best=None):\n",
    "        if tree is None:\n",
    "            return best\n",
    "        k = self.k\n",
    "        # 更新 branch\n",
    "        if point[0][depth % k] < tree.location[depth % k]:\n",
    "            next_branch = tree.left_child\n",
    "        else:\n",
    "            next_branch = tree.right_child\n",
    "        if not next_branch is None:\n",
    "            best = next_branch.location\n",
    "        return self._search(point,\n",
    "                            tree=next_branch,\n",
    "                            depth=depth + 1,\n",
    "                            best=best)\n",
    "\n",
    "    def fit(self, X):\n",
    "        self.kdtree = self._fit(X)\n",
    "        return self.kdtree\n",
    "\n",
    "    def predict(self, X):\n",
    "        res = self._search(X, self.kdtree)\n",
    "        return res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x点的最近邻点是(2, 3)\n"
     ]
    }
   ],
   "source": [
    "KNN = KdTree()\n",
    "X_train = np.array([[2, 3], [5, 4], [9, 6], [4, 7], [8, 1], [7, 2]])\n",
    "KNN.fit(X_train)\n",
    "X_new = np.array([[3, 4.5]])\n",
    "res = KNN.predict(X_new)\n",
    "\n",
    "x1 = res[0]\n",
    "x2 = res[1]\n",
    "\n",
    "print(\"x点的最近邻点是({0}, {1})\".format(x1, x2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "----\n",
    "参考代码：https://github.com/wzyonggege/statistical-learning-method\n",
    "\n",
    "本文代码更新地址：https://github.com/fengdu78/lihang-code\n",
    "\n",
    "习题解答：https://github.com/datawhalechina/statistical-learning-method-solutions-manual\n",
    "\n",
    "中文注释制作：机器学习初学者公众号：ID:ai-start-com\n",
    "\n",
    "配置环境：python 3.5+\n",
    "\n",
    "代码全部测试通过。\n",
    "![gongzhong](../gongzhong.jpg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
